Random spherical hyperbolic diffusion

TL;DR

This paper develops a series solution for hyperbolic diffusion on the sphere, analyzes its properties, and applies it to cosmic microwave background data, providing insights into the stochastic modeling of such phenomena.

Contribution

It derives an exact series solution for hyperbolic diffusion on the sphere with random initial conditions, including convergence analysis and applications to cosmic microwave background data.

Findings

Series expansion converges with quantifiable error bounds

Sample H"older continuity relates to angular power spectrum decay

Numerical results support theoretical convergence and applicability

Abstract

The paper starts by giving a motivation for this research and justifying the considered stochastic diffusion models for cosmic microwave background radiation studies. Then it derives the exact solution in terms of a series expansion to a hyperbolic diffusion equation on the unit sphere. The Cauchy problem with random initial conditions is studied. All assumptions are stated in terms of the angular power spectrum of the initial conditions. An approximation to the solution is given and analysed by finitely truncating the series expansion. The upper bounds for the convergence rates of the approximation errors are derived. Smoothness properties of the solution and its approximation are investigated. It is demonstrated that the sample H\"older continuity of these spherical fields is related to the decay of the angular power spectrum. Numerical studies of approximations to the solution and applications to cosmic microwave background data are presented to illustrate the theoretical results.

Equations

$$\mathbf{q}(\mathbf{x},t+{\tau }_0)\approx {\tau }_0\frac{\partial \mathbf{q}(\mathbf{x},\mathbf{t})}{\partial t}+\mathbf{q}(\mathbf{x},\mathbf{t})=-k\nabla u(\mathbf{x},t),$$

$$\frac{1}{c^2}\frac{{\partial }^{2}u(\mathbf{x},t)}{\partial t^2}+\frac{1}{D}\frac{\partial u(\mathbf{x},t)}{\partial t}={\nabla }^{2}u(\mathbf{x},t),$$

$$\frac{1}{c^2}\frac{{\partial }^2\phi (\mathbf{x},t)}{\partial t^2}-{\nabla }^2\phi (\mathbf{x},t)+\frac{m^2{c}^2}{{\hslash }^2}\phi (\mathbf{x},t)=0,$$

$$d{s}^2$$

$$g^{\mu \nu }{\nabla }_{\mu }{\nabla }_{\nu }\phi (\mathbf{x},t)+\left(\frac{m^2{c}^2}{{\hslash }^2}+\xi R\right)\phi (\mathbf{x},t)=0,$$

$$|\det g{|}^{-1/2}\frac{\partial }{\partial x^{\mu }}\left[|\det g{|}^{1/2}{g}^{\mu \nu }\frac{\partial \phi (\mathbf{x},t)}{\partial x^{\nu }}\right]+\left(\frac{m^2{c}^2}{{\hslash }^2}+\xi R\right)\phi (\mathbf{x},t)=0.$$

$$\frac{1}{c^2}\left[\frac{2a^{\prime }(\eta )}{a}\frac{\partial \phi (\mathbf{x},t)}{\partial \eta }+\frac{{\partial }^2\phi (\mathbf{x},t)}{\partial {\eta }^2}\right]-\sum _{i=1}^3\frac{{\partial }^2\phi (\mathbf{x},t)}{\partial x^i\partial x^i}=0.$$

$$\eta ={\eta }_1+\frac{R_0}{c_1}\ln \left(1+\frac{c_1[t-t_1]}{a_1{R}_0}\right),$$

$$a=a_1[1+e^{c_1(\eta -{\eta }_1)/R_0}],$$

$$\bar{\phi }(R,\theta ,\varphi ,\eta )=\frac{3}{R^3}{\int }_0^R{r}^2\phi (r,\theta ,\varphi ,\eta )dr,$$

$$\frac{1}{c^2}\left[\frac{2a^{\prime }(\eta )}{a}\frac{\partial }{\partial \eta }+\frac{{\partial }^2}{\partial {\eta }^2}\right]\bar{\phi }(R,\theta ,\varphi ,\eta )-\{\frac{1}{R^2\sin (\theta )}\frac{\partial }{\partial \theta }\left[\sin (\theta )\frac{\partial }{\partial \theta }\right]$$

$$+\frac{1}{R^2{\sin }^2(\theta )}\frac{{\partial }^2}{\partial {\varphi }^2}\}\bar{\phi }(R,\theta ,\varphi ,\eta )=\frac{3}{R}\frac{\partial \phi (R,\theta ,\varphi ,\eta )}{\partial R}$$

$$+\frac{2}{R^3}{\int }_0^R\left(\frac{1}{\sin (\theta )}\frac{\partial }{\partial \theta }\left[\sin (\theta )\frac{\partial \bar{\phi }(r,\theta ,\varphi ,\eta )}{\partial \theta }\right]+\frac{1}{{\sin }^2(\theta )}\frac{{\partial }^2\bar{\phi }(r,\theta ,\varphi ,\eta )}{\partial {\varphi }^2}\right)dr$$

$$-\underset{r\to 0}{\lim }{\left(\frac{r}{R}\right)}^3\left\{\frac{1}{r^2\sin (\theta )}\frac{\partial }{\partial \theta }\left[\sin (\theta )\frac{\partial \bar{\phi }(r,\theta ,\varphi ,\eta )}{\partial \theta }\right]+\frac{1}{r^2{\sin }^2(\theta )}\frac{{\partial }^2\bar{\phi }(r,\theta ,\varphi ,\eta )}{\partial {\varphi }^2}\right\}.$$

$$\frac{1}{c^2}\frac{{\partial }^2\bar{\phi }(R,\theta ,\varphi ,\eta )}{\partial {\eta }^2}+\frac{c_1}{c^2{R}_0}\frac{\partial \bar{\phi }(R,\theta ,\varphi ,\eta )}{\partial \eta }=\frac{1}{R_0^2}{\Delta }_{{\mathbb{S}}^2}\bar{\phi }(R,\theta ,\varphi ,\eta ).$$

$$\frac{{\partial }^2\bar{\phi }(R,\theta ,\varphi ,\hat{\eta })}{\partial {\hat{\eta }}^2}+\frac{\partial \bar{\phi }(R,\theta ,\varphi ,\hat{\eta })}{\partial \hat{\eta }}=\frac{1}{{\hat{R}}^2}{\Delta }_{{\mathbb{S}}^2}\bar{\phi }(R,\theta ,\varphi ,\hat{\eta }),$$

$$\hat{\eta }=\eta /{\eta }_s=\eta c^2/D=\eta c_1/R_0;\hat{R}=R_0/{\ell }_s=R_{0}c/D=c_1/c.$$

$$H^2=H_0^2({\Omega }_{\Lambda }+{\Omega }_{m}a(t{)}^{-3}+{\Omega }_{\gamma \nu }a(t{)}^{-4}+{\Omega }_{\kappa }a(t{)}^{-2}).$$

$$H_e=(2{\Omega }_m{)}^{1/2}{a}_e^{-3/2}{H}_0,c_e={\dot{a}}_e{a}_e{R}_0,D=\frac{c^2{R}_e}{c_e}=\frac{c^2}{H_e{a}_e}.$$

$${\mathbb{S}}^2=\left\{\mathbf{x}\in {\mathbb{R}}^3:\Vert \mathbf{x}\Vert =1\right\}\subset {\mathbb{R}}^3$$

$$\tilde{\sigma }(d\mathbf{x})=\sigma (d\theta ,d\phi )=\sin \theta d\theta d\phi ,(\theta ,\phi )\in [0,\pi )\times [0,2\pi ).$$

$$T=\left\{T(\theta ,\phi )=T_{\omega }(\theta ,\phi ):0\le \theta <\pi ,0\le \phi <2\pi ,\omega \in \Omega \right\},$$

$$T(\theta ,\phi )=\sum _{l=0}^{\infty }\sum _{m=-l}^l{a}_{lm}{Y}_{lm}(\theta ,\phi ),$$

$$\underset{L\to \infty }{\lim }\mathbf{E}\left(\underset{{\mathbb{S}}^2}{\int }{\left(T(\theta ,\phi )-\sum _{l=0}^L\sum _{m=-l}^l{Y}_{lm}(\theta ,\phi )a_{lm}\right)}^2\sin \theta d\theta d\phi \right)=0.$$

$$\underset{L\to \infty }{\lim }\mathbf{E}{\left(\tilde{T}(\mathbf{x})-\sum _{l=0}^L\sum _{m=-l}^l{\tilde{Y}}_{lm}(\mathbf{x})a_{lm}\right)}^2=0,$$

$${\tilde{Y}}_{lm}(\mathbf{x})=Y_{lm}(\theta ,\phi )=d_{lm}\exp (im\phi )P_l^m(\cos \theta ),$$

$$d_{lm}=(-1{)}^m{\left[\frac{(2l+1)(l-m)!}{4\pi (l+m)!}\right]}^{1/2},$$

$$P_l^m(x)=(-1{)}^m(1-x^2{)}^{m/2}\frac{d^m}{d{x}^m}{P}_l(x),P_l(x)=\frac{1}{2^{l}l!}\frac{d^l}{d{x}^l}(x^2-1{)}^l.$$

$${\int }_0^{\pi }{\int }_0^{2\pi }{Y}_{lm}^{\ast }(\theta ,\phi )Y_{l^{\prime }{m}^{\prime }}(\theta ,\phi )\sin \theta d\phi d\theta ={\delta }_l^{l^{\prime }}{\delta }_m^{m^{\prime }},$$

$$Y_{lm}^{\ast }(\theta ,\phi )=(-1{)}^m{Y}_{l(-m)}(\theta ,\phi ),$$

Full text

1 Introduction

The linear telegraph equation was introduced in the 1880s in Heaviside’s model of transmission lines (see, e.g. Berg ). Since then, the same linear partial differential equation has arisen in several quite different contexts. In the 1950s, Cattaneo Cattaneo introduced the hyperbolic heat equation that has a bounded speed of propagation of temperature disturbances, unlike the classical parabolic heat equation that has an unbounded propagation speed. A bounded speed of propagation also prohibits the generation of new cosmic structures that are correlated over space-like separated regions. The large-scale coherent structures that are observed in the cosmic microwave background, are understood to be remnants of acoustic waves in the plasma universe, seeded by a very short inflationary period of superluminal expansion (e.g. Dodelson , Weinberg ). Restricting the post-inflation propagation of disturbances to sub-luminal speeds is guaranteed by choosing the parameter $c$ in the hyperbolic heat equation to be less than or equal to the speed of light. This has some rationale in relativistic geometry, and it maintains the second law of thermodynamics Ali .

In Section 2, we add another application of the same equation. When $c$ is chosen to be the speed of light, the hyperbolic heat equation is indeed equivalent, by choosing an appropriate material coordinate system with conformal time coordinate, to the general relativistically covariant scalar Klein-Gordon equation, minimally coupled to an expanding space-time. In that coordinate system, the material radius of the expanding universe is constant.

Our focus in later sections is on the hyperbolic diffusion of random disturbances on the surface or interior of a sphere. In all of the aforementioned applications, this may be viewed as a canonical initial value problem on a compact manifold, as expressed in the classic texts such as Carslaw . However, unlike the classic texts, we are interested more in random initial conditions that include structures that cannot be causally connected. The evolution of fields and their correlations, under speed-limited diffusion, is of primary interest. In the case of the cosmic microwave temperature, typically with relative fluctuations of the order of $10^{-4}$, the currently available data are indeed represented on a spherical surface, with little reference to a radial coordinate, see PLANCK2 , PLANCK1 ,MarinucciPeccati11 .

Recent years have witnessed an enormous amount of attention, in the astrophysical and cosmological literature, on investigating spherical random fields. The empirical motivation for these studies comes from the current cosmological research. The NASA satellite mission WMAP and the ESA mission Planck, see PLANCK2 , PLANCK1 , probe Cosmic Microwave Background radiation (CMB) to an unprecedented accuracy. Figure 2 shows measurements of the CMB temperature intensity from Planck 2015 results used as an illustration in this paper. Figure 2 plots the corresponding best-fit scaled angular power spectrum. CMB can be viewed as a signature of the distribution of matter and radiation in the very early universe, and as such it is expected to yield very tight constraints on physical models for the Big Bang and subsequent phase transitions and nucleosynthesis. For the density fluctuations of this field, the highly popular inflationary scenario predicts a Gaussian distribution, whereas alternative cosmological theories, such as topological defects or non-standard inflationary models, predict otherwise.

The seminal work of Applegate et. al. Applegate showed that neutron-rich regions due to the longer mean free path of neutrons compared to protons, could strongly influence the distribution of helium and deuterium during nucleosynthesis. When the temperature cooled to allow weak nuclear re-combination and freeze-out, ($\frac{3}{2}kT<800$ keV), the neutron diffusion length was around 0.08 pc or 0.3 light-years Kurki . Some basic considerations on diffusion length are given in Appendix A. After allowing for neutron back-diffusion to proton-rich regions where neutrons have been further depleted by fusion Terasawa , the observed distribution of elements and temperature variations constrains the parameters of heterogeneous models Barrow , Iocco .

Fluctuations in CMB observations may also have a non-physical origin, i.e. they might be generated by systematic errors in the CMB map, such as noise which has not been properly removed, contamination from the galaxy or distortions in the optics of the telescope. A proper understanding of the density distributions of fluctuations is also instrumental for correct inference on the physical constants which can be estimated from CMB radiation.

From a mathematical point of view, properties of random fields defined by SPDE on Euclidean spaces is a well studied area, see, for example, Anh0 and the references therein. However the known results in the literature are not directly translatable to manifolds. Therefore, the problems of approximating and studying sample regularity of random fields on the sphere have attracted much recent attention, see, for example, Anh , Lang , MarinucciPeccati11 , Xiao and the references therein. It was shown that the convergence rate of approximation schemes based on truncated series expansions of such fields is often related to the decay of the angular power spectrum. Another line of investigations started in the paper Kozachenko . It studied approximations of Gaussian isotropic random fields on the sphere, but used different models and types of convergence. In this research we continue these investigations for spherical random fields generated by stochastic hyperbolic diffusion equations.

The discussion in the paper and the obtained results give some indication that the considered stochastic hyperbolic diffusion equation on the sphere can be used to capture important statistical and spectral properties of the CMB. The current research uses a different approach and a stochastic model comparing to Anh ,Kozachenko , Lang , and Xiao . As the approach yields the explicit solution in terms of series of elementary functions it could be very useful for various qualitative and numerical studies. The numerical studies suggest that the proposed approximations to solutions have an optimal order of convergence.

The paper is organised as follows. Section 2 and Appendix A provide physical motivation and detailed justification for stochastic models studied in the paper. Basic results and definitions about spherical isotropic random fields and their spectral and covariance representations are given in Section 3. Section 4 first derives the solution of the non-random hyperbolic diffusion equation on the sphere. Then, it gives the main results about solutions of the hyperbolic diffusion equation with random initial conditions. The convergence rates of truncated solutions to the exact solution are derived in Section 5. Smoothness properties of the solution and its approximation are studied in the next section. Section 7 and Appendix C present numerical studies and applications to CMB data. Short conclusions and future directions are discussed in Section 8. All proofs are given in Appendix B.

All numerical computations and simulations in this paper were performed using the software R version 3.5.3 and Python version 3.6.7. The results were derived using the HEALPix representation of the CMB data, see gor and http://healpix.sourceforge.net. In particular, the R package rcosmo Fryer 2 , Fryer was used for accessing the CMB data, computations and visualisations of the obtained results. The Python package healpy was used for fast spherical harmonics transformations of CMB maps. 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/

2 Derivation of hyperbolic diffusion equations

This section presents physical motivation and detailed justification for mathematical and stochastic models considered in the following sections. It also discusses and provides values of the model parameters for numerical studies in Section 7.

Conduction of heat results from energy transfer during collisions of constituent particles or, in rigid lattices, the transfer of photons of vibrational energy from one lattice site to another. These transfer processes have a typical time delay $\tau_{0}$, leading to a relaxation time $\tau_{0}$ in a continuum dynamical model. For example, particles in a Boltzmann gas of hard spheres of radius $r_{0}$ have mean kinetic energy $\frac{1}{2}mv^{2}=\frac{3}{2}k_{B}u$ and mean free path is $\ell=k_{B}u/4\sqrt{2}\pi r_{0}^{2}p$ where $k_{B}$ is Boltzmann’s constant, $u({\bf x},t)$ is absolute temperature, $v({\bf x},t)^{2}$ is mean squared velocity and $p({\bf x},t)$ is pressure Thompson . This means that the average time before a given particle collides is $\tau_{0}=\ell/v$. Before this time, the temperature gradient has no effect on the energy flux. This would lead to a delayed PDE but as demonstrated by Cattaneo Cattaneo , it is convenient to approximate a delayed differential equation by a higher-order differential equation that is local in time. From the lag between instantaneous heat flux density ${\bf q}({\bf x},t)$ and the gradient of instantaneous temperature $u({\bf x},t)$, there follows a higher-order correction to Fourier’s law simply by Taylor expansion:

[TABLE]

where ${\bf q}({\bf x},t),{\bf x}\in\mathbb{R}^{3},t>0,$ is a vector field, $u({\bf x},t),{\bf x}\in\mathbb{R}^{3},t>0,$ is a scalar field, $\tau_{0}$ and $k$ are constants.

This is a two-parameter generalisation of Fourier’s law, wherein $k$ is the thermal conductivity. Then by the local conservation of heat energy, $\rho C\frac{\partial u({\bf x},t)}{\partial t}=-\nabla\cdot{\bf q}({\bf x},t)$, where $\rho$ is density of the medium and $C$ is specific heat capacity. This leads directly to Cattaneo’s hyperbolic heat/diffusion equation,

[TABLE]

where $D=k/C\rho$ and $c=\sqrt{D/\tau_{0}},$ which is the least upper bound for the speed of propagation. This can easily be seen by constructing attenuated travelling sinusoidal wave solutions by separation of variables.

The simplest relativistic extension of the non-relativistic Schrödinger matter wave equation for a particle of mass $m$ is the Klein-Gordon equation in flat Minkowski space,

[TABLE]

which describes a spin-zero matter field Wigner such as a Higgs scalar Boson Higgs . In standard physics notation, $\hbar=h/2\pi$, where $h$ is Planck’s constant. The Klein-Gordon equation agrees with the Schrödinger equation at low energies (e.g. Bjorken ). The scalar field $\varphi=\varphi({\bf x},t),{\bf x}\in\mathbb{R}^{3},t>0,$ has an interpretation in quantum field theory after second quantisation when it is regarded as an operator rather than a classical function Schweber .

Now consider a scalar field minimally coupled to a spatially flat de Sitter universe, with expansion factor $a(t)$, depending on cosmic time $t$. The space-time metric is

[TABLE]

Here, $x^{i}$ are material coordinates. For example, these would remain constant at the centre of mass of a typically moving galaxy, where $t$ is the proper time. $\eta$ is the conformal time coordinate, $\eta=\int[1/a(t)]dt$, the adoption of which renders the metric to be conformal to the Minkowski metric. In this coordinate system $g_{\mu\nu}$ is a diagonal matrix with entries $g_{00}=c^{2}a^{2}$ and $g_{jj}=-a^{2}$ for $j=1,2,3.$ When adapted to the more general space-time, the scalar covariant Klein-Gordon wave equation is (e.g. Birrell )

[TABLE]

where $\nabla_{\mu}$ is the covariant derivative with respect to $x^{\mu}$ and $g^{\mu\nu}$ is the matrix inverse of $g_{\mu\nu}$.

We consider the case of negligible quantum particle mass $m$, and zero coupling ($\xi=0$) to the Ricci scalar $R$, in (3). This case is still called ”minimal coupling” because the space-time metric influences the covariant derivative and the Laplace-Beltrami operator. In PBPZ , it was convenient to analyse the unattenuated field $a\varphi({\bf x},t).$ In the current context, we retain as the field of interest, $\varphi({\bf x},t)$ which is attenuated as the universe expands.

It is well known (e.g. Birrell ) that the Laplace-Beltrami operator, which is the second-order operator acting in the first term of (3), can be expressed more conveniently in terms of the determinant of the metric tensor, $\det g=\det[g_{\mu\nu}]$. Then (3) is equivalent to

[TABLE]

For the expanding spherically symmetric universe with $\det g=-c^{2}a^{8}$, a direct calculation gives

[TABLE]

Near some time $t=t_{1}$ in the past, the expansion factor may be approximated by a linear function $a=a_{1}+c_{1}(t-t_{1})/R_{0}$, where $R_{0}$ is the current radius and $c_{1}$ is the speed of expansion at time $t_{1}$. In that case,

[TABLE]

[TABLE]

and $a^{\prime}(\eta)/a=\frac{c_{1}}{2R_{0}}$ at $t=t_{1}$. The covariant Klein-Gordon equation is approximated by the relativistic diffusion equation, which is the same as the hyperbolic diffusion equation (2) with $D=c^{2}R_{0}/c_{1}$. Up to the present time, CMB data is represented on a 2-sphere, with little reference to radial variation. This prompts one to take a volume average along the radial direction in a spherical sector with fixed small differential increment of solid angle, throughout (5). Defining

[TABLE]

where $(r,\theta,\phi)$ are material polar coordinates, (5) implies

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The two limits on the right hand side sum to $-\lim_{r\to 0}\left(\frac{r}{R}\right)^{3}\frac{1}{r^{2}}\Delta_{\mathbb{S}^{2}}\varphi,$ where $r^{-2}\Delta_{\mathbb{S}^{2}}\varphi$ is the surface Laplacian on a sphere of vanishingly small radius. Since the Laplacian is not expected to be more singular than $\mathcal{O}(r^{-2})$, this limit is zero. The free boundary condition of zero mass-energy-flux across the expanding maximal radius $a(\eta)R$, reduces in material coordinates to $\frac{\partial\varphi}{\partial r}=0$ at $r=R$. This annuls the first term on the right hand side of the above. The remaining integral term is $2/3$ times the sector radially doubly averaged divergence of $\nabla\varphi$. The sector average is weighted towards values at large $r$. Double averages are weighted even more towards values at large $r$. Doubly averaged fluctuations in $\varphi$ are expected to be small compared to singly averaged fluctuations. As an indicative study of a scalar field, we set the right hand side of (2) to zero. Then near $\eta=\eta_{1}$, the equation for $\bar{\varphi}$ reduces to

[TABLE]

Note that $R_{0}^{-2}\Delta_{\mathbb{S}^{2}}$ is the Laplace-Beltrami operator on a spherical surface of radius $R_{0}$. In the material coordinate system, the material radius is constant $R_{0}$ at all times. There is a unique rescaling of variables that further reduces (6) to normalised form,

[TABLE]

[TABLE]

The single remaining parameter $\hat{R}$ involves the expansion rate at time $t_{1}$ which is related to the Hubble parameter. In particular, we are interested in a time $t_{1}$ when the cosmic electromagnetic radiation began to propagate through a transparent recombined atomic medium, whose properties may still be inferred.

The cosmic background electromagnetic radiation became free for unimpaired propagation into the distant future, around the narrow band of time of atomic recombination. Seemingly by coincidence, this narrow period was close to the time of equality of matter density and radiation density expressed in the same units (e.g. Hir ). Henceforth, subscript zero will denote a parameter value measured at the current time in history whereas a subscript $e$ will denote its value at the time of equality. At the present time, the Hubble constant is measured to be $H_{0}=\dot{a}(t)/a(t)=0.7$ km/s/Mpc and the Hubble time is $1/H_{0}=1.4\times 10^{10}$ yr. From the Friedman equation, the Hubble ‘constant’ $H$ actually depends on time according to

[TABLE]

The various components $\Omega_{J}$ are fractions of current total energy that are dark energy $\Omega_{\Lambda}$, mass density $\Omega_{m}$, radiation energy density $\Omega_{\gamma\nu}$ and curvature energy density $\Omega_{\kappa}$. Within experimental error, the universe is spatially flat (e.g. NASA WMAP website NASA ), so the total mass-energy is the critical energy for the universe to be spatially flat. The present estimates, with $a(t)=1$, have dark energy $\Omega_{\Lambda}=0.714$, and matter $\Omega_{m}$=0.286 consisting of dark matter (0.240) and baryons (0.046). The component $\Omega_{\kappa}$ due to spatial curvature is measured to be zero and the component due to combined photon and neutrino radiation ($\Omega_{\gamma\nu}$) is negligible at the current time. However at other times, the components of mass-energy, relative to the current critical mass-energy, are multiplied by various powers of $a(t)$. For example around the age of 370,000 years after the big bang, within the short period of atomic recombination, the radiation and matter components $\Omega_{m}a(t)^{-3}$ and $\Omega_{\gamma\nu}a(t)^{-4}$ were equal and the other components were negligible. This occurs at expansion factor $a_{e}=\Omega_{\gamma\nu}/\Omega_{m}=4.15\times 10^{-5}$. Consequently,

[TABLE]

The conformal time scale appearing above is $\eta_{s}=D/c^{2}=R_{e}/c_{e}=1/H_{e}$, which is exactly the Hubble time at $t=t_{e}$. This is approximately $9\times 10^{-8}/H_{0}$, where $1/H_{0}$ is the current Hubble time. The material length scale evaluates to $\ell_{s}=\frac{4.15^{3/2}}{2\times 5.72^{1/2}}\times 10^{-7}c/H_{0}.$ This is approximately 0.02 times the radius of the universe at recombination time.

3 Isotropic random fields

This section introduces basic notations and background by reviewing some results in the theory of spherical random fields from the monograph Yadrenko (see, also Ivanov , Lang , Leonenko99 , MarinucciPeccati11 ).

We will use the symbol $C$ to denote constants which are not important for our exposition. Moreover, the same symbol may be used for different constants appearing in the same proof.

Consider a sphere in the three-dimensional Euclidean space

[TABLE]

with the Lebesgue measure (the area element on the sphere)

[TABLE]

A spherical random field on a complete probability space $(\Omega,\mathcal{F},\mathbf{P})$, denoted by

[TABLE]

in the spherical coordinate system, or $\widetilde{T}=\{\widetilde{T}(\mathbf{x})$ , $\mathbf{x}\in\mathbb{S}^{2}\}$ in the Cartesian coordinates, is a stochastic function defined on the sphere $\mathbb{S}^{2}.$

We consider a real-valued second-order spherical random field $T$ that is continuous in the mean-square sense. Note that MarinucciPeccati13 proved that the covariance function of a measurable finite-variance isotropic random field on the sphere is necessarily everywhere continuous.

Under these conditions, the field $T$ can be expanded in the mean-square sense as a Laplace series (see, Yadrenko , p. 73, Leonenko99 , p. 33, or MarinucciPeccati13 , p.123):

[TABLE]

where $\{Y_{lm}(\theta,\varphi)\}$ represents the complex spherical harmonics. The spectral representation (7) can be viewed as a Karhunen-Loève expansion, which converges in the Hilbert space $L_{2}(\Omega\times\mathbb{S}^{2},\sin\theta d\theta d\varphi)$, that is,

[TABLE]

According to the Peter-Weyl theorem (see MarinucciPeccati13 , p.69), the expansion (7) also converges in the Hilbert space $L_{2}(\Omega),$ for every $\mathbf{x}\in\mathbb{S}^{2}$, that is, for each $\mathbf{x}\in\mathbb{S}^{2},$

[TABLE]

where $\{\tilde{Y}_{lm}(\mathbf{x})\}$ represents the complex spherical harmonics of the Cartesian variable $\mathbf{x}.$

Recall that for $-l\leq m\leq l$ it holds

[TABLE]

[TABLE]

where $P_{l}^{m}(\cdot)$ denotes the associated Legendre polynomial with the indices $l$ and $m,$ and $P_{l}(\cdot)$ is the $l$-th Legendre polynomial, i.e.

[TABLE]

The spherical harmonics have the following properties

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\delta_{l}^{l^{\prime}}$ is the Kronecker delta function, the symbol * means the complex conjugation and $\mathbf{0}$ corresponds to $\varphi=\theta=0.$ The random coefficients $a_{lm}$ in the Laplace series (7) can be obtained through inversion arguments in the form of mean-square stochastic integrals

[TABLE]

As $T$ is real-valued, then, by the property (9), it holds

[TABLE]

The field $\widetilde{T}(\mathbf{x})$ is called isotropic (in the weak sense) on the sphere $\mathbb{S}^{2}$ if $\mathrm{E}\widetilde{T}(\mathbf{x})^{2}<\infty$ and its first and second-order moments are invariant with respect to the group $SO(3)$ of rotations in $\mathbb{R}^{3},$ i.e.

[TABLE]

for every $g\in SO(3)$ and $\mathbf{x},\mathbf{y}\in\mathbb{S}^{2}.$ This is equivalent to saying that the mean $\mathbf{E}T(\theta,\varphi)=c=constant$ (without loss of generality we assume that $c=0),$ and that the covariance function $\mathbf{E}T(\theta,\varphi)T(\theta^{\prime},\varphi^{\prime})$ depends only on the angular distance $\Theta=\Theta_{PQ}$ between the points $P=(\theta,\varphi)$ and $Q=(\theta^{\prime},\varphi^{\prime})$ on $\mathbb{S}^{2}.$

The field is isotropic if and only if

[TABLE]

Thus, $\mathbf{E}|a_{lm}|^{2}=C_{l},$ $m=0,\pm 1,...,\pm l.$ The series $\left\{C_{1},C_{2},...,C_{l},...\right\}$ is called the angular power spectrum of the isotropic random field $T(\theta,\varphi).$

From (7) and (11) we deduce that the covariance function of an isotropic random fields has the following representation

[TABLE]

where

[TABLE]

If $T(\theta,\varphi)$ is an isotropic Gaussian field, then the coefficients $a_{lm},$ $m=-l,\dots,l,$ $l\geq 1,$ are complex-valued independent Gaussian random variables if $m\not=-m^{\prime},$ with

[TABLE]

if $C_{l}>0$, or degenerate to zero if $C_{l}=0.$

4 Solution for stochastic spherical hyperbolic diffusion

First this section derives solutions of non-random hyperbolic diffusion equations. Then the obtained results are used in to obtain solutions of diffusion equations with random initial conditions.

Consider the following hyperbolic diffusion equation, also known as the telegraph equation (see Kolesnik ) or relativistic diffusion equation on sphere

[TABLE]

with the initial conditions

[TABLE]

where $\mathbf{x}=(x_{1},x_{2},x_{3})\in\mathbb{S}^{2},$ $c>0,\;D>0,$ and $k$ are some constants, and $\Delta_{\mathbb{S}^{2}}$ is the Laplace operator on the sphere $\mathbb{S}^{2}$ and $\delta(\mathbf{x})$ is the Dirac delta-function.

Note, that in the unit spherical coordinates equation (15) takes the form

[TABLE]

[TABLE]

where

[TABLE]

is the Laplace-Beltrami operator on the sphere.

It is known (see, i.e., MarinucciPeccati13 , p.72) that the eigenvalue problem for Laplace operator on the sphere has the following exact solution

[TABLE]

where $\{\tilde{Y}_{lm}(\mathbf{x})\}$ is the system of spherical harmonics. Therefore, it is natural to seek a solution of the problem (15)-(16) in the form of the series

[TABLE]

where

[TABLE]

$\tilde{\sigma}(d\mathbf{x})=\sin\theta d\theta d\varphi$.

The proof of the following result is given in Appendix B.

Theorem 4.1

The solution $\tilde{p}(\mathbf{x},t)$ to the hyperbolic diffusion point-source initial value problem (15)-(16) is given by the series

[TABLE]

[TABLE]

where $\mathbf{1}_{\{\cdot\}}$ denotes the binary indicator function,

[TABLE]

[TABLE]

Now we use the results of Theorem 4.1 to derive solutions of diffusion equations with random initial conditions. The random field $u(\theta,\varphi,t)$ is defined by the following hyperbolic diffusion equation on the sphere

[TABLE]

[TABLE]

where $\Delta_{(\theta,\varphi)}$ is the Laplace-Beltrami operator on the sphere given by (17).

Now, the random initial conditions are determined by the Gaussian isotropic random field on the sphere

[TABLE]

[TABLE]

where $a_{lm},$ $m=-l,\dots,l,$ $l\geq 0,$ are complex-valued independent Gaussian random variables satisfying (12) and (14).

Theorem 4.2

If the angular power spectrum {$C_{l},l=0,1,2,...\}$ of the random field $T(\theta,\varphi)$ from the initial condition (23) satisfies assumption (13), then the random solution $u(\theta,\varphi,t)$ of the initial value problem (22)-(24) is given by the convergent in the Hilbert space $L_{2}(\Omega\times\mathbb{S}^{2},\sin\theta d\theta d\varphi)$ random series

[TABLE]

where

[TABLE]

are stochastic processes with

[TABLE]

and

[TABLE]

Moreover, its covariance function is given by

[TABLE]

[TABLE]

The proof of Theorem 4.2 is given in Appendix B.

5 Convergence study of approximate solutions

The results in Section 4 provide a series representation of the random field $u(\theta,\varphi,t).$ To investigate contributions of different terms in the representation it is important to study their finite cumulative sums. In this section it is done by analysing truncated series expansions of the solution $u(\theta,\varphi,t)$ of the initial value problem (22)-(24). We demonstrate the role of the decay rate of the angular power spectrum. These results are also important for various numerical studies. In particular, they can be used to determine the required number of terms in the truncated series to get a specified accuracy of the approximate solutions.

The approximation $u_{L}(\theta,\varphi,t)$ of truncation degree $L\in\mathbb{N}$ to the solution $u(\theta,\varphi,t)$ in (25) is defined by

[TABLE]

for $\theta\in[0,\pi),$ $\varphi\in[0,2\pi),$ $t>0.$

The following theorem gives the convergence rate of the approximation $u_{L}(\theta,\varphi,t)$ to the solution $u(\theta,\varphi,t).$ The order of the convergence rate is determined by the high frequency magnitudes of the angular spectrum.

Theorem 5.1

Let $u(\theta,\varphi,t)$ be the solution to the initial value problem (22)-(24) and $u_{L}(\theta,\varphi,t)$ be the approximation of truncation degree $L\in\mathbb{N}$ of $u(\theta,\varphi,t).$

Then, for $t>0$ the truncation error is bounded by

[TABLE]

Moreover, for $L>\frac{\sqrt{D^{2}k^{2}+c^{2}}-Dk}{2Dk}$ it holds

[TABLE]

where the constant $C$ depends only on the parameters $c$, $D$ and $k.$

The proof of Theorem 5.1 is given in Appendix B and does not depend on $\theta$ and $\varphi.$ Therefore, the statement also holds in $L_{2}(\Omega)$ norm over the sphere ${\mathbb{S}}^{2}.$ Hence, we obtain the following results.

Corollary 1

Uniformly over $\theta\in[0,\pi),$ $\varphi\in[0,2\pi)$ the results of Theorem 5.1 are also valid for the mean squared truncation error

[TABLE]

[TABLE]

Corollary 2

Let the angular power spectrum $\{C_{l},l=0,1,2,...\}$ of the random field $T(\theta,\varphi)$ from the initial condition (23) decay algebraically with order $\alpha>2,$ that is, there exist constants $C>0$ and $l_{0}\in\mathbb{N}$ such that $C_{l}\leq C\cdot l^{-\alpha}$ for all $l\geq l_{0}.$

Then, the approximation $u_{L}(\theta,\varphi,t)$ converges to the solution $u(\theta,\varphi,t)$ of the initial value problem (22)-(24) and

1. (i)

for $L>\max(l_{0},\frac{\sqrt{D^{2}k^{2}+c^{2}}-Dk}{2Dk})$ the truncation error is bounded by

[TABLE] 2. (ii)

for any $\varepsilon>0$ it holds

[TABLE] 3. (iii)

for all $\theta\in[0,\pi),$ $\varphi\in[0,2\pi)$ and $t>0$ the truncation error is asymptotically almost surely bounded by

[TABLE]

where $\beta\in\left(0,\frac{\alpha-3}{2}\right)$ and $\alpha>3.$

The proof of Corollary 2 is given in Appendix B.

6 Hölder continuity of solutions and their truncated approximations

In this sections we investigate properties of the solution $u(\theta,\varphi,t)$ of the initial value problem (22)-(24) and its approximation $u_{L}(\theta,\varphi,t).$ We show that they are Hölder continuous fields. It is demonstrated how the decay of the angular power spectrum is related to the Hölder continuity in mean square of the corresponding random field.

It follows from MarinucciPeccati13 that the field $u(\theta,\varphi,t)$ is mean square continuous. However, obtaining sample Hölder continuity of this field requires stronger assumptions on the decay of the angular power spectrum of $T(\theta,\varphi)$ than (13).

Theorem 6.1

Let $u(\theta,\varphi,t)$ be the solution to the initial value problem (22)-(24) and the angular power spectrum {$C_{l},l=0,1,2,...\}$ of the random field $T(\theta,\varphi)$ from the initial condition (23) satisfies the assumption

[TABLE]

Then there exists a constant C such that for all $t>0$ it holds

[TABLE]

where the constant $C$ depends only on the parameters $c$, $D$ and $k.$

Replacing $\sum_{l=0}^{\infty}$ by $\sum_{l=L}^{\infty}$ in the proof of Theorem 6.1 in Appendix B we obtain the Hölder continuity of the approximations to the solution.

Corollary 3

Let the assumptions of Theorem 6.1 hold true. Then there exists a constant $C_{L}$ such that for all $t>0$ it holds

[TABLE]

where the constant $C_{L}$ depends only on the parameters $c$, $D$ and $k.$

The following result is proven in Appendix B. It provides an upper bound on $p$th moments of the solution increments in time.

Corollary 4

Let the assumptions of Theorem 6.1 hold true. Then, for each $p>0,$ there exists a constant C such that for all $t>0$ it holds

[TABLE]

where the constant $C$ depends only on the parameters $p,$ $c$, $D$ and $k.$

Finally, we present continuity properties of the solution at time $t$ with respect to the geodesic distance on the sphere.

Corollary 5

Let $u(\theta,\varphi,t)$ be the solution to the initial value problem (22)-(24) and the angular power spectrum {$C_{l},l=0,1,2,...\}$ of the random field $T(\theta,\varphi)$ from the initial condition (23) satisfies the assumption

[TABLE]

Then, there exists a constant C such that for all $t>0$ it holds

[TABLE]

where $\Theta$ is the angular distance between $(\theta,\varphi)$ and $(\theta^{\prime},\varphi^{\prime})$ and the constant $C$ depends only on the parameters $c$, $D$ and $k.$

7 Numerical studies

In this section, we present detailed numerical studies of the solution $u(\theta,\varphi,t)$ of the initial value problem (22)-(24) and its approximation $u_{L}(\theta,\varphi,t).$ We investigate the convergence rates of the approximation to the solution and evolutions of the solution, the covariance function and its power spectrum over time.

We use the data with measurements of the CMB temperature intensity from Planck 2015 results, see PLANCK2 and PLANCK1 . Figure 2 shows the CMB map produced from the SMICA (a component separation method for CMB data processing) pipeline data at $N_{\rm side}=1024$ at $10$ arcmin resolution with $12,582,912$ HEALPix points. Figure 2 plots the best-fit LCDM scaled angular power spectrum $D_{l}=l(l+1)C_{l}/(2\pi),$ $l=2,...,2508$, of the CMB map at the recombination time. The scaled CMB angular power $D_{l}$ is shown as a function of the harmonic number $l.$ It begins at $l=2$ as for $l=0$ and 1 it can not be reliably estimated using only $2l+1$ values.

We use the coefficients $a_{lm}$ and the angular power spectrum of CMB temperature intensities in Figures 2 and 2 as the initial condition of the Cauchy problem (22)-(24). For numerical studies and R computations in this section we assume that the angular spectrum of the random field $T(\theta,\varphi)$ is vanished if $l$ is greater than 2508. Thus, we use $u_{L_{0}}(\theta,\varphi,t)$ with $L_{0}=2508$ as a substitution of the solution $u(\theta,\varphi,t).$ Also, as the temperature of the ensemble of decoupled photons has continued to diminish and now shows very small variability in the range $2.7260\pm 0.0013$ K, in plots we use the same colour scheme but different scales compared to the intensity map in Figure 2 that corresponds to time $t=0$ in the model (22)-(24). It helps better visualise the solutions and their approximation errors.

7.1 Evolution of solutions

Study of the evolution of CMB field is critical to unveil important properties of the present and primordial universe PLANCK2 , Dodelson . The CMB map can be modelled as a realization of the random field $u(\theta,\varphi,t)$ on $\mathbb{S}^{2}.$ We demonstrate its evolution due to the model (22)-(24). For the following numerical examples we will use the random field $T(\theta,\varphi)$ with the angular power spectrum of CMB temperature intensities given in Figure 2.

In the case of the model (22)-(24) it follows from (25) and (26) that the angular spectrum evolution over time is determined by the multiplication factor $\exp\left(-{c^{2}t}/{2D}\right)[A_{l}^{2}(t)+B_{l}^{2}(t)].$ Since the attenuation factor is $\exp(-c^{2}t/2D)$ the dimensionless time $t^{\prime}=c^{2}t/2D$ will be used in all following plots. Figure 4 shows the original scaled angular power spectrum $D_{l}$ in red and the angular power spectra at time $t^{\prime}=0.02$ and 0.04 in green and blue respectively. It is observed that $D_{l}$ changes little from the original values over short periods. The deviations increase with increasing $l$ which is consistent with the changes of the multiplication factor shown in Figure 4 and cosmological theories showing that higher multipoles are changing faster.

Figure 6 displays the covariance $\mathbf{Cov}(u(\theta,\varphi,0),u(\theta^{\prime},\varphi^{\prime},t^{\prime}))$ defined by (29). Three lines depict covariances at the time lags $t^{\prime}=0,$ 0.02, and 0.04 as functions of the angular distance $\Theta.$ To further understand the impact of time and the angular distance on the covariance we produce 3d-plots showing the covariance as a function of the time lag $t^{\prime}$ and the angular distance $\Theta,$ see Figure 6. The plots in Figures 6 and 6 are normalised by dividing each value by the variance at time 0 and the angular distance 0, i.e. by $\mathbf{Cov}(u(0,0,0),u(0,0,0)).$ It is observed that the covariance decays very rapidly. It changes very little over short time periods, except $\Theta$ values close to 0. As the angular power spectrum decreases very quickly only its values at small multipoles have the principal impact on covariances. Relative large changes of $D_{l}$ at high frequencies (see Figure 4) do not substantially change the covariance function. It is evidenced from Figure 6 where the three lines almost coincide. Hence, it would not be reasonable to use the covariance function to characterise fine changes in CMB maps over short periods of time. However, as the correlations change little over time it can help in studying the CMB at earlier epochs.

We investigated realisations of the solution that correspond to scaled angular power spectra for $c=1,$ $D=1$ and $k=0.01$ at early conformal times. Their temperature maps were rather similar to the original map in Figure 2 that corresponds to time $t^{\prime}=0.$ The temperature field becomes smoother and its range narrows when time increases. Figure 8 depicts differences of the temperature intensities at time $t^{\prime}=0$ and 0.04. The fine structure of changes is shown in Figure 8 that visualises a small triangular ”equatorial” region around the Milky Way galaxy. No regions with extreme changes compare to other locations can be seen in the both figures.

Figure 10 shows a realisation of the solution that corresponds to the case $c=1,$ $D=1$ and $k=0.05$ at time $t^{\prime}=0.04.$ The temperature intensity differences between the solution fields with $k=0.01$ and $k=0.05$ (Figure 10) at time $t^{\prime}=0.04$ are visualised in Figure 10. As expected the higher value 0.05 of the diffusivity parameter $k$ results in a ”blurred” realisation of the map with $k=0.01.$ The difference field does not exhibit any specific spatial pattern.

7.2 Convergence rates of approximations to solutions

Now we analyse approximations to solutions and corresponding approximation errors depending on the truncation degree $L.$

We analysed the spatial error fields $u(\theta,\varphi,t^{\prime})-u_{L}(\theta,\varphi,t^{\prime})$ of approximations $u_{200}(\theta,\varphi,0.04)$ and $u_{400}(\theta,\varphi,0.04)$ to the solution $u(\theta,\varphi,0.04)$ with $c=1,$ $D=1$ and $k=0.01$ at time $t^{\prime}=0.04.$ The approximation error field for the case $L=200$ was rather similar to the true map which indicates that more terms are required to reconstruct fine details of the temperature intensity fields. For the case $L=400$ the error field did not exhibit any specific spatial pattern. An increase in the approximation accuracy was also evidenced by the decrease of the mean squared error from $5.988295\cdot 10^{-09}$ to $4.075033\cdot 10^{-09}.$

Figure 11 shows the difference of the mean ${L_{2}(\Omega\times{\mathbb{S}}^{2})}$ truncation errors \Big{\|}u(\theta,\varphi,t^{\prime})-u_{L}(\theta,\varphi,t^{\prime})\Big{\|}_{L_{2}(\Omega\times{\mathbb{S}}^{2})} and their upper bounds (30) in Theorem 5.1 on a natural logarithmic scale. The case of the SPDE (22) with $c=1,$ $D=1$ and $k=0.1$ at $t^{\prime}=10$ is considered. The difference is plotted as a function of $L.$ The plot confirms that both the error and its upper bound asymptotically vanish when $L$ increases. Moreover, the convergence rates of the error and its upper bound are of the same order and differ only by a constant multiplication factor. Thus, the approximations to the solutions achieve an optimal order of convergence. Figure 11 also suggests that the convergence rate is faster than the power one.

Some additional numerical studies of the dependence on time and the parameters are given in Appendix C. The code in the folder ”Research materials” from https://sites.google.com/site/olenkoandriy/ can be used to obtain maps of realisations for different parameters that were not included in this section or to experiment with other combinations of the parameters.

8 Conclusions and future work

This research provides physical motivation and justification for stochastic diffusion models in CMB studies. The Cauchy problem with random initial conditions for hyperbolic diffusion equations on the unit sphere was considered. Properties of exact and approximate solutions were investigated. The numerical studies illustrated the obtained results and examined the sensitivity of solutions to parameters using the CMB data. They indicated that the model is flexible enough to capture some of the statistical properties of the CMB.

It was shown that properties of the solutions are determined by the decay of the angular power spectrum. The analysis of the CMB covariance revealed that dependencies between CMB observations rapidly decrease with angular distance between their locations. The numerical studies demonstrated that the CMB evolution under this SPDE model results in most significant changes of CMB temperature spectrum at the high frequency range in the first place. At the same time relatively large changes of $D_{l}$ at high CMB frequencies have little impact on the covariance functions and CMB maps that remain almost unchanged. The contribution of high frequency components to the CMB field is decreasing with $l$ at a faster rate than a power-law.

The sensitivity analysis to parameters at the lower frequency spectrum range demonstrated that changes in time and the parameters $c$ and $D$ have a substantial impact on the error, while $k$ showed almost no effect.

It would be interesting to extend the obtained results and

• (1)

apply this model to backward studies of CMB;

• (2)

compare this model with other possible evolution scenarios;

• (3)

investigate 3-dimensional SPDEs as the next generation CMB experiment, CMB-S4, will be collecting 3D spatial data;

• (4)

explore applications to other spherical data in physics and earth sciences.

Acknowledgements.

This research was supported under the Australian Research Council’s Discovery Project DP160101366. N. Leonenko was supported in part by Cardiff Incoming Visiting Fellowship Scheme, International Collaboration Seedcorn Fund, Data Innovation URI Seedcorn Fund. We are also grateful for the use of data of the Planck/ESA mission from the Planck Legacy Archive. The authors are also grateful to the referees for their careful reading of the paper and suggestions that helped to improve the paper.

Appendix A: Diffusion length of a local disturbance.

Consider a density disturbance $u$ of total mass $Q$ originating at the origin. The well-known point source solution to linear diffusion in three space dimensions is given by

[TABLE]

The density level set at some low significance value $u$ is at

[TABLE]

so

[TABLE]

The level set reaches its maximum extent when $\frac{dr}{dt}=0$, implying $t=\frac{(Q/u)^{2/3}}{8e\pi D},$ so the diffusion length is

[TABLE]

For example, for a mass disturbance the size of a solar mass, and a neutron diffusion length of 0.3 light-year at the temperature of neutrino dissociation from weak nuclear interactions, estimated from Applegate and Kurki , the marginal disturbance density $u$ is around 1 solar mass per cubic light year. This is meant to have occurred at a time when the cosmological expansion factor $a(t)$ was less than $10^{-3}$, so after expansion to the current level, the equivalent marginal density would be less than one nucleon mass per cubic metre, around the current mean density of the universe.

Appendix B: Proofs.

Proof of Theorem 4.1. By substituting (19) into equation (15) and using (18), we obtain

[TABLE]

To find particular solutions of (31), we need to solve the ordinary differential equation

[TABLE]

The initial conditions for this equation can be determined from (20) and (16) and they are

[TABLE]

The characteristic equation of (32) is $\frac{1}{c^{2}}z^{2}+\frac{1}{D}z+l(l+1)k^{2}=0,$ with the roots $z_{1,2}=-{c^{2}}/(2D)\pm K_{l}.$ Therefore, the general solution of equation (32) is given by the formula:

[TABLE]

where $M_{1},M_{2}$ are some constants. From the initial conditions in (33) we obtain

[TABLE]

Thus, the solution of the Cauchy problem (32)-(33) is given by

[TABLE]

[TABLE]

Returning now to (19), we obtain the solution of the Cauchy problem (15)-(16) in the form

[TABLE]

[TABLE]

Note that the multiplier of $Q_{l}(\mathbf{x})$ on the right-hand side of (34) equals

[TABLE]

By substituting this expression into (34), we get

[TABLE]

Finally, using $K_{l}^{\prime}$ and rewriting the Green function we obtain the statement of the theorem.

Proof of Theorem 4.2. The solution of the initial value problem (22) - (24) can be written as a spherical convolution of the Green function $p(\theta,\varphi,t)$ from Section 4 and the random field $T(\theta,\varphi),$ if the corresponding Laplace series converges in the Hilbert space $L_{2}(\Omega\times\mathbb{S}^{2},\sin\theta d\theta d\varphi).$

Let the two functions $f_{1}(\cdot)$ and $f_{2}(\cdot)$ on the sphere $\mathbb{S}^{2}$ belong to the space $L_{2}(\mathbb{S}^{2},\sin\theta d\theta d\varphi)$ and have the Fourier-Laplace coefficients

[TABLE]

Recall (see, i.e., Dunkel ) that their non-commutative spherical convolution is defined as the Laplace series

[TABLE]

with the Fourier-Laplace coefficients given by

[TABLE]

provided that the series (35) converges in the corresponding Hilbert space.

Thus, the random solution $u(\theta,\varphi,t)$ of equation (22) with the initial values determined by (23) and (24) can be written as a spherical random field with the following Laplace series representation

[TABLE]

provided that this series is convergent in the Hilbert space $L_{2}(\Omega\times\mathbb{S}^{2},\sin\theta d\theta d\varphi),$ where $p_{t}=p(\theta,\varphi,t)$ is given by Theorem 4.1 and $T$ is given by (23). The complex Gaussian random variables $a_{lm}^{(t)}$ are given by

[TABLE]

where $a_{l0}^{(p_{t})}=Y_{l0}^{\ast}(\mathbf{0})d_{l}(\theta,\varphi,t)$ and

[TABLE]

[TABLE]

It gives the first statement of the theorem.

By the addition formula for spherical harmonics (see, i.e., MarinucciPeccati13 , p.66)

[TABLE]

where $P_{l}(\cdot)$ is the $l$-th Legendre polynomial (see (8)), and $\cos\Theta$ is the angular distance between the points $(\theta,\varphi)$ and $(\theta^{\prime},\varphi^{\prime})$ on $\mathbb{S}^{2}.$

Using (10) we obtain that the random field $u(\theta,\varphi,t)$ is isotropic if and only if the covariance structure of the solution (25) can be written in the form

[TABLE]

[TABLE]

which gives the result in (29) provided the series (29) converges for every fixed $t$ and $t^{\prime},$ that is

[TABLE]

Noting that $\left|P_{l}(\cos\Theta)\right|\leq 1,$ only a finite number of terms $A_{l}$ is non-zero, and there is a constant $C$ such that $\sup_{t\geq 0}|B(t)|<C,$ we obtain that condition (38) follows from (13). This condition on the angular spectrum $C_{l},l\geq 0,$ guarantees the convergence of the series (36) in the Hilbert space $L_{2}(\Omega\times{\mathbb{S}}^{2},\sin\theta d\theta d\varphi).$

Proof of Theorem 5.1. The approximation $u_{L}(\theta,\varphi,t)$ is a centered Gaussian random field, i.e. $\mathbf{E}u_{L}(\theta,\varphi,t)=0$ for all $L\in\mathbb{N},$ $\theta\in[0,\pi),$ $\varphi\in[0,2\pi),$ and $t>0.$ Therefore,

[TABLE]

[TABLE]

[TABLE]

By (27) and (28) we get

[TABLE]

Hence, for all $L\in\mathbb{N}$ it holds

[TABLE]

For $l>\frac{\sqrt{D^{2}k^{2}+c^{2}}-Dk}{2Dk}$ it follows from (27) that $A_{l}(t)\equiv 0$. Therefore, by (39) and (40) we obtain

[TABLE]

Proof of Corollary 2. The statement (i) immediately follows from (30) and the estimate

[TABLE]

Then, applying Chebyshev’s inequality, we get the upper bound in (ii).

Finally, (iii) follows from statement (ii) and the Borel–Cantelli lemma as

[TABLE]

Proof of Theorem 6.1. Let $h$ belong to a bounded neighbourhood of the origin. It follows from (14), (26), (27), (28) and (37) that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We start by showing how to estimate the first summand in (41). By (27), for the case $l=0$ we obtain

[TABLE]

[TABLE]

For $l>0$ we will use the upper bound

[TABLE]

[TABLE]

[TABLE]

By properties of $\cosh(\cdot)$ and $\sinh(\cdot)$ we get

[TABLE]

Then, applying (27) and noting that only a finite number of $A_{l}$ is non-vanished (namely, only if $l\in\left[0,\frac{\sqrt{D^{2}k^{2}+c^{2}}-Dk}{2Dk}\right]$) we obtain the following estimates

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Now we estimate the second summand in (41) as

[TABLE]

[TABLE]

[TABLE]

Using (40) and applying the inequalities $|\cos(x)-\cos(y)|\leq 2\left|\sin\left(\frac{x-y}{2}\right)\right|\leq|x-y|$ and $|\sin(x)-\sin(y)|\leq|x-y|$ we obtain

[TABLE]

[TABLE]

Note that for all $l\geq 0$ it holds

[TABLE]

Applying the above estimates to (41) we obtain

[TABLE]

[TABLE]

which completes the proof.

Proof of Corollary 4. Note that $u(\theta,\varphi,t)$ is a centered Gaussian random field and for any centered Gaussian random variable $X$ it holds

[TABLE]

Applying this result to the statement of Theorem 6.1 we obtain

[TABLE]

Proof of Corollary 5. By(29) it holds

[TABLE]

[TABLE]

[TABLE]

Applying the next property of Legendre polynomials (see, for example, Lang , p.16) $|1-P_{l}(x)|\leq 2|1-x|^{\gamma}(l(l+1))^{\gamma},\quad\gamma\in[0,1],$ and the upper bounds (40), we obtain that uniformly in $t\geq 0$

[TABLE]

Appendix C: Sensitivity to parameters.

To further understand the impact of time and the model parameters on the difference of the mean ${L_{2}(\Omega\times{\mathbb{S}}^{2})}$-errors and their upper bound (30) we produced 3d-plots showing the difference as a function of the truncation degree $L$ and each parameter provided that other parameters are fixed. These plots are displayed in Figures 15, 15, 15, and 15.

In all cases the difference between the error and its upper bound asymptotically vanish when $L$ increases. Figure 15 demonstrates that the difference is a decreasing function of time $t^{\prime},$ which is expected as the series representation (25) of the solutions $u(\theta,\varphi,t^{\prime})$ has the multiplication factor $\exp\left(-t^{\prime}\right)=\exp\left(-{c^{2}t}/(2D)\right)$ exponentially decaying in time. The differences are extreme at the origin and decrease when time or the parameter $c$ increases, see Figure 15. For the parameter $D$ the situation depicted in Figure 15 is opposite and the difference is increasing in $D$ which is expected as the multiplication factor is exponentially decaying in $D^{-1}.$ Finally, Figure 15 suggests that the parameter $k$ seems have no substantial impact on the difference.