1 Introduction

The NASA mission WMAP and the European Space Agency mission Planck have collected highly accurate cosmological data, resulting in a precise map depicting the distribution of cosmic microwave background radiation (CMB) [2, 3]. It is expected that new experiments, in particular, within CMB-S4 Collaboration and the Euclid mission, will provide CMB measurements at new unprecedented precision. The necessity of modelling and analysing them have recently attracted increasing attention to spherical random fields. From the mathematical point of view and for modelling purposes, CMB maps can be regarded as a single realization of a random field on a sphere of a large radius. The sphere plays a role of the underlying space and expands in time. For the dark energy-dominated era, the expansion factor has the exponential form [8]. This expansion impacts the stochastic diffusion and should be incorporated in a model for the evolution of CMB. Contrary to the other models in the literature, this paper takes into consideration this expansion of space–time, which leads to a stochastic diffusion equation with non-constant coefficients that evolve with the expansion factor.

The CMB spectrum indicates that since the last scattering around 380,000 years after the big bang, the universe has been transparent to electromagnetic radiation. However the universe is not transparent to charged cosmic ray particles. They are deviated by magnetic fields as they pass close to galaxies. Recent estimates of the number of galaxies in the observable universe range from \(2\times 10^{12}\) to \(6\times 10^{12}\). In any angular aperture of observation, there will be part of at least one galaxy. Extragalactic cosmic ray particles reach the earth in a small number of showers each year. They are distinguished from local galactic cosmic rays by their high particle energy, greater than \(5 \times 10^{18} eV\). This corresponds to charged particles typically arriving at more than half the speed of light. They have likely been travelling for a very long time in cosmic terms, during which they would have been deviated by a number of galaxies. This results in a long-term diffusive redistribution of matter throughout the universe. This process is very complex. Effective diffusion occurs partly by scattering due to magnetic fields [18] and also by motion of the magnetic field lines themselves [22]. It has also been argued that the temperature gradients formed from galaxies can effectively repel particles, enough to avoid gravitational capture [33].

We refer to the monograph [27] for the systematic exposition of the main results of the theory of spherical random fields and their cosmological applications. The paper [23] studied isotropic random fields on high-dimensional spheres and established connections between the smoothness of the covariance kernel approximation rates and the decay of their angular power spectrum.

One important direction of the modern theory of random fields is the exploration of extremes of random fields including fields given on manifolds; see the classical results in [13, 31, 35] and their modern generalizations. The monograph [11] examined extremes of sub-Gaussian fields while the expected Euler characteristic method developed by Adler and Taylor was demonstrated in [4]. The publication [12] provided the asymptotics of excursion probabilities for both smooth and non-smooth spherical Gaussian fields. Excursion probabilities for spherical sub-Gaussian random fields were studied in [32]. Another approach utilized limit theorems for sojourn measures, see [25] and [26].

Stochastic partial differential equations (SPDE) are the main tool to model the evolution of spherical random fields over time. They have been extensively studied, see [15, 23, 30] and the references therein. Several models were recently presented in [5, 9, 10, 24] that employed stochastic hyperbolic diffusion equations and modelled various types of evolution of spherical random fields.

This article integrates the aforementioned approaches and extends them to the context of SPDEs within an expanding space–time. The equations studied in this paper differ from the mentioned SPDEs and are given as a hyperbolic diffusion on the expanding sphere. Motivated by cosmological applications, the exponential expansion is used in the considered model, which leads to a stochastic diffusion equation with non-constant coefficients. The article also conducts the numerical analysis of the solution of the considered model. The numerical analysis section investigates the evolution of the solution of the studied SPDE over space–time, the structure of its space–time dependencies and its extremes. The CMB intensity map from the mission Planck and its spectrum are used as initial conditions. The numerical analysis confirms and visualizes the obtained theoretical results.

The main novelties of the paper include:

  • The consideration of diffusion within an expanding space–time framework;

  • The examination of equations with non-constant coefficients dependent on the expansion factor values;

  • An exploration of both the local and asymptotic properties of the solutions and their respective approximations;

  • An analysis of excursion probabilities associated with the solutions and their approximations.

The paper is structured as follows: Sect. 2 provides the main definitions and notations. Section 3 presents the diffusion model within an expanding time–space universe. Then, this section investigates the initial value problem for this diffusion equation and explores the properties of the derived non-random solutions. Section 3 is dedicated to the examination of the equation with random initial conditions. The solution to the equation and its associated covariance function are derived. Section 5 investigates properties of stochastic solution and its corresponding approximations. Section 6 studies excursion probabilities related to the solution and its approximations. Finally, Sect. 7 presents numerical study that illustrate the properties of the solution and the obtained results.

2 Main definitions and notations

This section reviews the main definitions and notations used in this paper and provides required background knowledge from the theory of random fields.

The notation \({\varvec{x}} = (x_1,x_2,x_{3})\in {\mathbb {S}}^2= \{\textbf{x}: \ ||{\textbf{x}}|| = 1,\ {\textbf{x}}\in {\mathbb {R}}^3\}\) is used to define Euclidean coordinates of points, while the notation \((\theta , \varphi ), \ 0\le \theta \le \pi , \ 0\le \varphi <2\pi ,\) is used for the corresponding spherical coordinates. In what follows, \( \Theta \) denotes the angular distance between points with spherical coordinates \((\theta , \varphi )\) and \((\theta ', \varphi ').\) \(||\cdot ||\) is the Euclidean norm in \({\mathbb {R}}^3,\) and \(Leb(\cdot )\) stands for the Lebesgue measure on the unit sphere \( {\mathbb {S}}^2.\) C with subindices represents a generic finite positive constant, the values of which are not necessarily the same in each appearance and may be changed depending on the expression.

By \(L^2({\mathbb {S}}^2)\), we denote the Hilbert space of square integrable functions on \({\mathbb {S}}^2\) with the following canonical inner product

$$\begin{aligned} \langle f,g \rangle _{L^2({\mathbb {S}}^2)} = \int \limits _0^\pi \int \limits _0^{2\pi } f(\theta ,\varphi )g^*(\theta ,\varphi )\sin \theta d\theta d \varphi ,\ \ \ f,g\in L^2({\mathbb {S}}^2), \end{aligned}$$

where \(g^*(\cdot )\) denotes a complex conjugate of \(g(\cdot )\) and the induced norm \( ||f||^2_{L^2({\mathbb {S}}^2)} = \langle f, f \rangle .\)

The spherical harmonics \(Y_{lm}(\theta ,\varphi ), \ l=0,1 \ldots , m = -l,..l,\) are given as

$$\begin{aligned} Y_{lm}(\theta ,\varphi ) = d_{lm}e^{ im\varphi }P_l^m(\cos \theta ), \end{aligned}$$

where

$$\begin{aligned}d_{lm} = (-1)^m\sqrt{\frac{2l+1}{4\pi }\frac{(l-m)!}{(l+m)!}},\end{aligned}$$

\(P_l^m(\cdot )\) are the associated Legendre polynomials with the indices l and m,  and \(P_l(\cdot )\) is the lth Legendre polynomial. The spherical harmonics \(Y_{lm}(\theta ,\varphi ), \ l=0,1, \ldots , m = -l,..l,\) form an orthogonal basis in the Hilbert space \(L^2({\mathbb {S}}^2).\)

The addition formula for the spherical harmonics states that

$$\begin{aligned} \sum _{m=-l}^lY_{lm}(\theta , \varphi ) Y^*_{lm}(\theta ', \varphi ') = \frac{2l+1}{4\pi }P_l(\cos \Theta ). \end{aligned}$$
(1)

By \(L^2(\Omega \times {\mathbb {S}}^2)\), we denote the Hilbert space of spherical random fields that have a finite norm

$$\begin{aligned} ||T||_{L^2(\Omega \times {\mathbb {S}}^2)}= E\left( \int \limits _{0}^{\pi }\int \limits _{0}^{2\pi } T^2(\theta ,\varphi ) \sin \theta d\varphi d\theta \right) .\end{aligned}$$

An isotropic spherical random field \(T(\theta ,\varphi ) \in L^2(\Omega \times {\mathbb {S}}^2)\) allows a representation as the following Laplace series [23]

$$\begin{aligned} T(\theta , \varphi ) = \sum _{l=0}^\infty \sum _{m=-l}^l a_{lm} Y_{lm}(\theta ,\varphi ), \end{aligned}$$
(2)

where the convergence is in the space \(L^2(\Omega \times {\mathbb {S}}^2).\) The random variables \(a_{lm}, \ l=0,1, \ldots ,\ m=-l, \ldots ,l,\) are given by the next stochastic integrals defined in the mean-square sense

$$\begin{aligned} a_{lm} = \int \limits _{0}^\pi \int \limits _{0}^{2\pi }T(\theta ,\varphi )Y^*(\theta ,\varphi )\sin \theta \textrm{d}\theta \textrm{d}\varphi . \end{aligned}$$
(3)

If \(T(\theta ,\varphi )\in L^2(\Omega \times {\mathbb {S}}^2)\) is a centred real-valued Gaussian random field, then it allows the representation (2), where \(a_{lm},\ l=0,1, \ldots , \ m=-l, \ldots ,l,\) are Gaussian random variables such that

$$\begin{aligned} a_{lm} = (-1)^ma_{l(-m)}, \quad Ea_{lm}=0, \quad Ea_{lm}a_{l'm'}^* = \delta _m^{m'}\delta _l^{l'}C_l. \end{aligned}$$

The sequence \(\{ C_l, l=0,1, \ldots \}\) is called the angular power spectrum of the isotropic random field \(T(\theta ,\varphi ).\) The series (2) converges in the \(L_2(\Omega \times {\mathbb {S}}^2)\) sense if it holds true [23, Section 2] that

$$\begin{aligned} \sum \limits _{l=0}^\infty C_l(2l+1)<+\infty . \end{aligned}$$
(4)

We assume that the condition (4) remains valid throughout all subsequent sections of the paper.

3 Spherical diffusion within expanding space–time

This section provides a justification for the model and obtains certain properties of solutions to the non-random version of the considered diffusion problem.

The universe is observed to have spatial cross sections with zero curvature. The Friedmann–Lemaître–Robertson–Walker metric (FLRW) on the sphere of radius r is

$$\begin{aligned} \textrm{d}s^2 = c^2 \textrm{d}t^2 - a^2(t)(r^2\textrm{d}\theta ^2 + r^2\sin ^2\theta \textrm{d}\varphi ^2),\end{aligned}$$

in the spherical space coordinates \((\theta , \varphi )\) and the cosmic time coordinate t,  where the term \(a(\cdot )\) is the expansion factor. Note that for the dark energy-dominated era [19], the expansion factor has the exponential form of the maximally symmetric de Sitter universe, \(a(t) = e^{\sqrt{\frac{\Lambda }{3}}ct}\) [8, 37].

The spherical diffusion is given by the following equation

$$\begin{aligned} \frac{1}{D} \frac{\partial {\widetilde{u}}}{\partial t} + g^{\mu \nu } \nabla _\mu \nabla _{\nu } {\widetilde{u}} =0, \end{aligned}$$
(5)

where \(\nabla _\mu \) is the covariant derivative operator and \(g^{\mu \nu }\) are the elements of the contravariant metric tensor, where the indices take values from 0 to 2. We also impose the following initial conditions

$$\begin{aligned} {\widetilde{u}}(t, \theta , \varphi )|_{t = 0} = \delta (\theta , \varphi ), \ \ \ \frac{\partial {\widetilde{u}}(t, \theta , \varphi )}{\partial t}{\biggl |}_{t = 0} =0. \end{aligned}$$
(6)

The covariant derivatives do not commute when they act on vectors, see [7]. However, the Laplace–Beltrami operator on the sphere can be expanded as follows

$$\begin{aligned} g^{\mu \nu }\nabla _\mu \nabla _\nu {\widetilde{u}} = \frac{1}{\sqrt{|g|}}\partial _\mu \sqrt{|g|}g^{\mu \nu } \partial _\nu {\widetilde{u}}, \end{aligned}$$
(7)

where g is the covariant metric tensor. For the above defined FLRW metric, g is the diagonal matrix with the elements \(c^{2}, -r^2a^{2}(t), -r^2a^{2}(t)sin^{2}\theta ,\) and the contravariant metric tensor is the inverse of g.

Changing to a conformal time coordinate \( \eta = \int \limits _0^t\frac{1}{a(s)}ds\) the FLRW metric becomes

$$\begin{aligned} \textrm{d}s^2 = e^2(\eta )(c^2 \textrm{d}\eta ^2 -r^2\textrm{d}\theta ^2 - r^2\sin ^2\theta \textrm{d}\varphi ^2),\end{aligned}$$

where \(e(\eta ):=\frac{\eta _\infty }{\eta _\infty -\eta }, \ \eta _\infty = \frac{1}{c}\sqrt{\frac{3}{\Lambda }},\) and the covariant metric tensor g becomes the diagonal matrix with the elements \(c^{2}e^2(\eta ), -r^2e^2(\eta ), -r^2e^2(\eta )\sin ^{2}\theta .\)

Using (7) and the above expression of the covariant metric tensor in terms of the conformal time, Eq. (5) can be written in the following form

$$\begin{aligned} \left( \frac{e(\eta )}{D} + \frac{e'}{c^2e}\right) \frac{\partial {\widetilde{u}}}{\partial \eta } + \frac{1}{c^2} \frac{\partial ^2 {\widetilde{u}}}{\partial \eta ^2} - \frac{1}{r^2} \frac{\partial ^2{\widetilde{u}}}{\partial \theta ^2} - \frac{1}{r^2 \sin ^2\theta }\frac{\partial ^2{\widetilde{u}}}{\partial \varphi ^2} -\frac{\cot \theta }{r^2}\frac{\partial {\widetilde{u}}}{\partial \theta } = 0. \end{aligned}$$
(8)

As for \(t = 0\) it holds that \(\eta = 0\) and \(\frac{\textrm{d}\eta }{\textrm{d}t}\big |_{t=0} = 1,\) the initial conditions (6) become

$$\begin{aligned} {\widetilde{u}}(\eta , \theta , \varphi )|_{\eta = 0} = \delta (\theta , \varphi ), \ \ \ \frac{\partial {\widetilde{u}}(\eta , \theta , \varphi )}{\partial \eta }{\biggl |}_{\eta = 0} =0. \end{aligned}$$
(9)

Equation (8) is different compared to the models studied in the literature, see [5, 9, 10, 24] and the references therein. Namely, for the underlying FLRW space–time metric the coefficient of the term \(\frac{\partial {\widetilde{u}}}{\partial \eta }\) in (8) is not constant and depends on the evolution of the expansion factor \(e(\eta ).\)

Theorem 1

The solution \({\widetilde{u}}(\eta , \theta , \varphi )\) of Eq. (8) with the initial conditions (9) is given by the following series

$$\begin{aligned}\sum _{l=0}^{\infty } F_l(\eta ) \sum _{m=-l}^l Y_{lm}^*({\textbf{0}}) Y_{lm}(\theta , \varphi ),\end{aligned}$$

where

$$\begin{aligned} F_l(\eta ) = {\left\{ \begin{array}{ll} 1, \ l =0,\\ (\eta _\infty - \eta )^\nu \big (K_1^{(l)}J_\nu (z_l(\eta _\infty - \eta ))+K_2^{(l)}Y_{\nu }(z_l(\eta _\infty - \eta ))\big ), \ l\in {\mathbb {N}}, \end{array}\right. }\ \end{aligned}$$
(10)
$$\begin{aligned} K_1^{(l)} = \frac{\pi z_l Y_{\nu -1}(z_l\eta _\infty )}{2 \eta _\infty ^{\nu -1}}, \ \ \ K_2^{(l)} = -\frac{\pi z_l J_{\nu -1}(z_l\eta _\infty )}{2 \eta _\infty ^{\nu -1}},\end{aligned}$$

\(z_l = \frac{c\sqrt{l(l+1)}}{r},\) \(\nu = \frac{c^2\eta _\infty }{2D}+1\) and 0 denotes a point on the unit sphere with the spherical coordinates \(\theta = 0, \ \varphi =0.\) \(J_\nu (\cdot )\) and \(Y_\nu (\cdot ),\) \(\nu \in {\mathbb {R}},\) are the Bessel functions of the first and second kind respectively.

Proof

Let \({\widetilde{u}} = L(\theta )Z(\varphi )F(\eta )\) be the solution of (8), then, after multiplying the equation by \(\frac{r^2sin^2\theta }{{\widetilde{u}}},\) one gets

$$\begin{aligned} r^2 \sin ^2\theta E(\eta ) \frac{F'}{F} + \frac{r^2 \sin ^2\theta }{c^2}\frac{F''}{F} - \sin ^2\theta \frac{L''}{L} - \sin \theta \cos \theta \frac{L'}{L} = \frac{Z''}{Z} = -m^2,\end{aligned}$$

where \(E(\eta ) = \frac{e(\eta )}{D} + \frac{e'}{c^2e},\) and m is a separation constant. Thus, \(Z(\varphi )=e^{im\varphi }\) and as it is \(2\pi \)-periodic, m is an integer.

Now separating the independent variable \(\theta \), one gets

$$\begin{aligned}r^2 E(\eta ) \frac{F'}{F} + \frac{r^2}{c^2}\frac{F''}{F}= \frac{L''}{L} + \textrm{cot}\theta \frac{L'}{L} - \frac{m^2}{\sin ^2\theta } = -l(l+1).\end{aligned}$$

By substituting \(x = cos\theta ,\) this equation for \(\theta \) is equivalent to the next associated Legendre equation [6, page 648]

$$\begin{aligned} (1-x^2)\frac{\textrm{d}^2L}{\textrm{d}x^2} -2x \frac{\textrm{d}L}{\textrm{d}x} -\frac{m^2 L}{1-x^2} +l(l+1)L = 0.\end{aligned}$$

The solution of the above equation has a singularity if \(l\notin {\mathbb {Z}}.\) If \(l=0,1, \ldots ,\) and \(m=-l, \ldots ,l,\) then the general solution is \(C_1 P_l^m(\textrm{cos} \theta ) + C_2 Q_l^m(\textrm{cos} \theta ),\) where \(P_l^m(\cdot )\) and \(Q_l^m(\cdot )\) are the associated Legendre polynomial and the Legendre function of the second kind. As \(Q_l^m(x)\) is singular at \(x = \pm 1,\) \(L(\theta ) = P_l^m(\cos \theta ).\) Due to the identity \(P^m_{-l} = P^m_{l-1},\ l=1,2, \ldots ,\) the solutions for negative integers l coincide with the obtained above solutions.

As \(e(\eta ) = \frac{\eta _\infty }{\eta _\infty -\eta },\) the equation for the separated independent variable \(\eta \) is

$$\begin{aligned} F'' + \frac{\eta _\infty c^2+D}{D(\eta _\infty - \eta )}F' + \frac{l(l+1)c^2}{r^2} F = 0. \end{aligned}$$
(11)

Let us denote the solutions to the above equation for \(l = 0,1, \ldots ,\) by \(F_l(\eta ).\) By straightforward calculations, one gets \(F_0 = K_1^{(0)}(\eta _\infty -\eta )^{2\nu } + K_2^{(0)},\ \nu = \frac{\eta _\infty c^2}{2D}+1,\) where the superscript (0) means that the constants \(K_1^{(0)}\) and \(K_2^{(0)}\) correspond to the solution of (11) with \(l = 0\).

By setting \(x=\eta _\infty -\eta ,\) one transforms (11) to

$$\begin{aligned} \frac{\textrm{d}^2F}{\textrm{d}x^2} + \left( -\left( \frac{\eta _\infty c^2}{D}+2\right) + 1\right) \frac{1}{x}\frac{\textrm{d}F}{\textrm{d}x} + \frac{l(l+1)c^2}{r^2} F = 0.\end{aligned}$$

Then, by [36, page 97] for \(l\in {\mathbb {N}}\) one gets

$$\begin{aligned} F_l = (\eta _\infty - \eta )^\nu \left( K_1^{(l)} J_\nu (z_l(\eta _\infty - \eta )) + K_2^{(l)} Y_{\nu }(z_l(\eta _\infty - \eta ))\right) .\end{aligned}$$

As \(Y_{lm}(\theta , \varphi ) = d_{lm}e^{im\phi }P_l^m(cos\theta ),\) the general solution of (8) is

$$\begin{aligned} \sum _{l=0}^\infty F_l(\eta ) \sum _{m = -l}^l d_{lm}^{-1}Y_{lm}(\theta , \varphi ), \end{aligned}$$
(12)

where

$$\begin{aligned} F_l(\eta ) = {\left\{ \begin{array}{ll} K_1^{(0)}(\eta _\infty - \eta )^{2\nu } + K_2^{(0)}, \ l =0,\\ (\eta _\infty - \eta )^\nu \big (K_1^{(l)}J_\nu (z_l(\eta _\infty - \eta ))+K_2^{(l)}Y_{\nu }(z_l(\eta _\infty - \eta ))\big ), \ l\in {\mathbb {N}}. \end{array}\right. }\ \end{aligned}$$
(13)

Due to the spherical harmonic closure relations [29, 1.17.25], the initial conditions (9) are equivalent to the following system of conditions

$$\begin{aligned} {\left\{ \begin{array}{ll} d_{lm}^{-1} F_l\big |_{\eta =0} = Y_{lm}^*({\textbf{0}}),\\ F_l'\big |_{\eta =0} = 0, \end{array}\right. }\ \end{aligned}$$
(14)

for all \(l = 0,1, \ldots \)

For \(l=0\), one easily gets that \(K_1^{(0)} = 0\) and \(K_2^{(0)} = d_{00}Y_{00}^*({\textbf{0}}).\) For \(l\in {\mathbb {N}}\), the conditions (14) become

$$\begin{aligned} {\left\{ \begin{array}{ll} (\eta _\infty - \eta )^{\nu }\big (K_1^{(l)} J_\nu (z_l(\eta _\infty - \eta )) + K_2^{(l)}Y_{\nu }(z_l(\eta _\infty - \eta ))\big )\big |_{\eta =0} = d_{lm}Y_{lm}^*({\textbf{0}}),\\ \left( (\eta _\infty - \eta )^{\nu }\big ( K_1^{(l)} J_\nu (z_l(\eta _\infty - \eta )) + K_2^{(l)}Y_{\nu }(z_l(\eta _\infty - \eta ))\big ) \right) '\big |_{\eta =0} = 0. \end{array}\right. }\ \end{aligned}$$

The last system is equivalent to

$$\begin{aligned} {\left\{ \begin{array}{ll} A_1^{(l)}K_1^{(l)} + B_1^{(l)}K_2^{(l)} = d_{lm}Y_{lm}^*({\textbf{0}}), \\ A_2^{(l)}K_1^{(l)} + B_2^{(l)}K_2^{(l)} = 0 \end{array}\right. }\ \end{aligned}$$

where

$$\begin{aligned} A_1^{(l)} = \eta _\infty ^\nu J_\nu (z_l\eta _\infty ), \ \ \ A_2^{(l)} = - \nu \eta _\infty ^{\nu -1}J_\nu (z_l\eta _\infty ) - \frac{z_l}{2}\eta _\infty ^\nu (J_{\nu -1}(z_l\eta _\infty ) - J_{\nu +1}(z_l\eta _\infty )),\\ B_1^{(l)} = \eta _\infty ^\nu Y_{\nu }(z_l\eta _\infty ), \ \ \ B_2^{(l)} = - \nu \eta _\infty ^{\nu -1}Y_{\nu }(z_l\eta _\infty ) - \frac{z_l}{2}\eta _\infty ^\nu (Y_{\nu -1}(z_l\eta _\infty ) - Y_{\nu +1}(z_l\eta _\infty )).\end{aligned}$$

Thus,

$$\begin{aligned} K_1^{(l)} = \frac{B_2^{(l)}d_{lm}Y_{lm}^*({\textbf{0}})}{A_1^{(l)}B_2^{(l)}-A_2^{(l)}B_1^{(l)}}, \ \ K_2^{(l)} = - \frac{A_2^{(l)}d_{lm}Y_{lm}^*({\textbf{0}})}{A_1^{(l)}B_2^{(l)}-A_2^{(l)}B_1^{(l)}}. \end{aligned}$$

After straightforward algebraic manipulations, one can see that

$$\begin{aligned} A_1^{(l)}B_2^{(l)}-A_2^{(l)}B_1^{(l)}= & {} \frac{z_l}{2}\eta _\infty ^{2\nu }\bigg ( -\big (J_{\nu +1}(z_l\eta _\infty )Y_{\nu }(z_l\eta _\infty ) - J_{\nu }(z_l\eta _\infty )Y_{\nu +1}(z_l\eta _\infty )\big )\\ {}{} & {} - \big (J_{\nu }(z_l\eta _\infty )Y_{\nu -1}(z_l\eta _\infty ) - J_{\nu -1}(z_l\eta _\infty )Y_{\nu }(z_l\eta _\infty ) \big )\bigg ).\end{aligned}$$

Using the Wronskian expression \(W(J_\nu (x), Y_{\nu }(x)) = J_{\nu +1}(x)Y_{\nu }(x) - J_{\nu }(x)Y_{\nu +1}(x) = \frac{2}{\pi x}\) (see [29, 10.5]), one gets

$$\begin{aligned} A_1^{(l)}B_2^{(l)}-A_2^{(l)}B_1^{(l)} = - \frac{2\eta _\infty ^{2\nu -1}}{\pi },\end{aligned}$$

from which follows that

$$\begin{aligned} K_1^{(l)} = \frac{\pi d_{lm} \bigg (\nu Y_{\nu }(z_l\eta _\infty )+\frac{z_l\eta _\infty }{2}(Y_{\nu -1}(z_l\eta _\infty )-Y_{\nu +1}(z_l\eta _\infty ))\bigg )Y_{lm}^*({\textbf{0}})}{2 \eta _\infty ^{\nu }}, \end{aligned}$$
$$\begin{aligned} K_2^{(l)} = -\frac{\pi d_{lm} \bigg (\nu J_\nu (z_l\eta _\infty )+\frac{z_l\eta _\infty }{2}(J_{\nu -1}(z_l\eta _\infty )-J_{\nu +1}(z_l\eta _\infty ))\bigg )Y_{lm}^*({\textbf{0}})}{2 \eta _\infty ^{\nu } }.\end{aligned}$$

By subsequently applying the identities \(Y_{\nu }'(x) = \frac{1}{2}(Y_{\nu -1}(x) - Y_{\nu +1}(x))\) and \(xY_\nu '(x) = xY_{\nu -1}(x)-\nu Y_{\nu }(x)\), one gets \(\nu Y_\nu (x) + \frac{x}{2}(Y_{\nu -1}(x)-Y_{\nu +1}(x)) = x Y_{\nu -1}(x).\) Thus,

$$\begin{aligned} K_1^{(l)} = \frac{\pi d_{lm} z_l Y_{\nu -1}(z_l\eta _\infty )Y_{lm}^*({\textbf{0}})}{2 \eta _\infty ^{\nu -1}}. \end{aligned}$$

Analogous transformations for the functions \(J_\nu (\cdot )\) lead to

$$\begin{aligned} K_2^{(l)} = -\frac{\pi d_{lm} z_l J_{\nu -1}(z_l\eta _\infty )Y_{lm}^*({\textbf{0}})}{2 \eta _\infty ^{\nu -1}}.\end{aligned}$$

By putting the above constants into (13) and (12), one finishes the proof of the theorem.\(\square \)

The following results derive some properties of the functions \(F_l(\eta ),\) which will be used later.

Lemma 1

For a fixed \(\eta \in [0,\eta _\infty )\) and any constant \(a>1\), it holds true that

$$\begin{aligned} Y_{a-1}(z_l\eta _\infty ) J_{a}(z_l(\eta _\infty -\eta )) - J_{a-1}(z_l\eta _\infty )Y_{a}(z_l(\eta _\infty -\eta )) = \frac{2\cos (z_l\eta )}{\pi z_l \sqrt{\eta _\infty (\eta _\infty -\eta )}} + O_\eta (z_l^{-2}), \end{aligned}$$
(15)

as \( l\rightarrow \infty ,\) where the terms \(O_\eta (z_l^{-2})\) may depend on \(\eta .\)

Proof

For the Bessel’s function of the first kind \(J_{a}(x)\) with \(a\ge -\frac{1}{2}\), the following holds true uniformly in \(x\ge 0\)

$$\begin{aligned} \left| J_a(x) - \left( \frac{2}{\pi x} \right) ^{\frac{1}{2}} \cos \left( x-\frac{1}{2}a \pi - \frac{1}{4}\pi \right) \right| \le \frac{d_a}{x^{\frac{3}{2}}}, \end{aligned}$$

where \(d_a\) is a constant depending on a,  see [28, Theorem 4.1].

An approximation for \(Y_\nu (x)\) analogously follows from the relationship [36, Section 3.6]

$$\begin{aligned} Y_a(x) = \frac{H_a^{(1)}(x) - H_a^{(2)}(x)}{2i}, \end{aligned}$$

where \(H_a^{(1)}(x)\) and \(H_a^{(2)}(x)\) are the Hankel functions of the first and second kind respectively. A simple modification of the proof of [28, Theorem 4.1] gives that, uniformly in \(x\ge 2,\) for the Bessel’s function of the second kind \(Y_{a}(x),\) \(a\ge -{1}/{2},\) it holds true

$$\begin{aligned} \left| Y_a(x) - \left( \frac{2}{\pi x} \right) ^{\frac{1}{2}}\sin \left( x-\frac{1}{2}a \pi - \frac{1}{4}\pi \right) \right| \le \frac{d_a}{x^{\frac{3}{2}}}. \end{aligned}$$
(16)

Thus, left-hand side of (15) is asymptotically equal to

$$\begin{aligned}{} & {} \frac{2}{\pi z_l \sqrt{\eta _\infty (\eta _\infty - \eta )}}\left( \sin \left( z_l \eta _\infty - \frac{(a-1)\pi }{2} - \frac{\pi }{4}\right) \cos \left( z_l(\eta _\infty -\eta ) - \frac{a\pi }{2} - \frac{\pi }{4}\right) \right. \\ {}{} & {} \quad \left. - \cos \left( z_l \eta _\infty - \frac{(a-1)\pi }{2} - \frac{\pi }{4}\right) \sin \left( z_l (\eta _\infty -\eta )- \frac{a\pi }{2} - \frac{\pi }{4}\right) \right) + O_\eta (z_l^{-2})\\ {}{} & {} \quad = \frac{2\cos (z_l\eta )}{\pi z_l \sqrt{\eta _\infty (\eta _\infty -\eta )}} + O_\eta (z_l^{-2}). \end{aligned}$$

\(\square \)

Lemma 2

For a fixed \(\eta \in [0,\eta _\infty )\), the following asymptotic behaviour holds true

$$\begin{aligned} F_l(\eta ) = \frac{\cos (z_l \eta )}{e^{\nu -1/2}(\eta )} + O_\eta (z_l^{-1}), \quad l\rightarrow \infty , \end{aligned}$$

where the terms \(O_\eta (z_l^{-1})\) may depend on \(\eta .\)

Proof

Let us consider separately the following expression

$$\begin{aligned}{} & {} K_1^{(l)}J_\nu (z_l(\eta _\infty - \eta ))+K_2^{(l)}Y_{\nu }(z_l(\eta _\infty - \eta )) \\ {}{} & {} \quad = \frac{\pi z_l }{2 \eta _\infty ^{\nu -1}}(Y_{\nu -1}(z_l\eta _\infty )J_\nu (z_l(\eta _\infty - \eta )) - J_{\nu -1}(z_l\eta _\infty ) Y_{\nu }(z_l(\eta _\infty - \eta ))).\end{aligned}$$

According to Lemma 1, the above expression equals to

$$\begin{aligned} \frac{\sqrt{e(\eta )}\cos (z_l\eta )}{\eta _\infty ^\nu } + O_\eta (z_l^{-1}). \end{aligned}$$

Thus,

$$\begin{aligned} F_l(\eta ) = \frac{(\eta _\infty - \eta )^\nu \sqrt{e(\eta )}\cos (z_l\eta )}{\eta _\infty ^\nu } + O_\eta (z_l^{-1}) =\frac{\cos (z_l \eta )}{e^{\nu -1/2}(\eta )} + O_\eta (z_l^{-1}). \end{aligned}$$

\(\square \)

Lemma 3

The functions \(F_l(\eta )\) are uniformly bounded on \(\eta \in [0,\eta _\infty )\) and \(\ l\in 0\cup {\mathbb {N}}.\)

Proof

Note that \(\nu >1.\) First, let us consider the absolute value of the first summand in (10)

$$\begin{aligned} |z_l(\eta _\infty -\eta )^\nu Y_{\nu -1}(z_l\eta _\infty ) J_\nu (z_l(\eta _\infty -\eta ))|. \end{aligned}$$

It is bounded by some constant C for all values of \(\eta \in [0,\eta _\infty )\) and \(\ l\in 0\cup {\mathbb {N}}.\) Indeed, \(|(z_l(\eta _\infty -\eta ))^\frac{1}{2} J_\nu (z_l(\eta _\infty -\eta ))|<C,\) as \(|\sqrt{x}J_\nu (x)|, x\in [0,\infty ),\) is a bounded function [28]. The terms \(|z_l^\frac{1}{2}Y_{\nu -1}(z_l\eta _\infty )|\) are uniformly bounded for all l as the function \(\sqrt{x} Y_\nu (x)\) is bounded on \([C,\infty ]\) for any fixed \(C>0.\) The later boundedness follows from the continuity of \(\sqrt{x} Y_\nu (x)\) on \([C,\infty ]\) and the inequality (16).

The absolute value of the second summand in (10) can be estimated as

$$\begin{aligned} |z_l (\eta _\infty -\eta )^\nu J_{\nu -1}(z_l\eta _\infty )Y_\nu (z_l(\eta _\infty -\eta ))|\le C |z_l^{\frac{1}{2}}(\eta _\infty -\eta )^\nu Y_\nu (z_l(\eta _\infty -\eta ))|, \end{aligned}$$
(17)

which is obtained by using the boundedness of the function \(|\sqrt{x}J_\nu (x)|, \ x\in [0,\infty ).\) The boundedness of the function on the right-hand side of (17) follows from the asymptotic behaviour \(Y_\nu (x)\sim \frac{C}{x^\nu },\) \(x\rightarrow 0,\) see [1, 9.1.11], (16) and the continuity of the function \(Y_\nu (x)\) on the interval \(x\in [C,\infty ).\) \(\square \)

4 Solution for stochastic spherical diffusion equation

In this section, we consider the problem from Sect. 3 with the initial conditions determined by an isotropic Gaussian random field.

Namely, we consider the following initial condition and its spherical harmonics expansion

$$\begin{aligned}{} & {} {u}(\eta , \theta , \varphi )|_{\eta = 0} = T(\theta ,\varphi ) = \sum _{l=0}^{\infty }\sum _{m=-l}^la_{lm}Y_{lm}(\theta ,\varphi ), \end{aligned}$$
(18)
$$\begin{aligned}{} & {} \frac{\partial {u}(\eta , \theta , \varphi )}{\partial \eta }{\biggl |}_{\eta = 0} =0, \end{aligned}$$
(19)

where \(a_{lm},\ l=0,1, \ldots , \ m=-l, \ldots ,l,\) are Gaussian random variables. Without loss of generality, we assume that the random field \(T(\theta , \varphi )\) is centred \(ET(\theta ,\varphi ) = 0.\)

Lemma 4

If the angular power spectrum \(\{C_l, l=0,1,\ldots \}\) of the Gaussian isotropic random field \(T(\theta , \varphi )\) satisfies the condition

$$\begin{aligned} \sum _{l=1}^\infty C_l l^{10}(2l+1) < +\infty , \end{aligned}$$
(20)

then \(T(\theta ,\varphi )\in C^2({\mathbb {S}}^2)\) a.s.

Proof

It follows from (1) that

$$\begin{aligned} Cov(T(\theta , \varphi ), T(\theta ', \varphi ')) = \sum _{l=0}^\infty C_l\sum _{m=-l}^lY_{lm}(\theta , \varphi ) Y^*_{lm}(\theta ', \varphi ')= (4\pi )^{-1} \sum _{l=0}^\infty C_l (2l+1)P_l(\cos \Theta ).\end{aligned}$$

Therefore, the statement of the lemma directly follows from [12, Lemma 3.3].\(\square \)

Lemma 4 provides the sufficient conditions for the random field \(T(\theta ,\varphi )\) to be a.s. twice continuously differentiable with respect to \(\theta ,\varphi .\) In the following results, we assume that \(T(\theta ,\varphi )\) is a.s. twice continuously differentiable or that the conditions of Lemma 4 hold true.

Lemma 5

Let \(K_T(\Theta )\) be a covariance function of an isotropic a.s. continuous Gaussian random field \(T(\theta ,\varphi ).\) There exists \(\Theta '>0\) such that \(K_T(\Theta )=K_T(0),\) for all \(\Theta \le \Theta ',\) if and only if \(T(\theta ,\varphi )=\xi , \ \xi \in L^2(\Omega ),\) for all \(\theta \) and \(\varphi ,\) with probability 1.

Proof

As the sufficiency is straightforward, let us proceed with the necessity.

Note that for any two spherical points \((\theta ,\varphi )\) and \((\theta ',\varphi ')\) at the angular distance \(\Theta ,\) it holds

$$\begin{aligned}E(T(\theta ,\varphi )- T(\theta ',\varphi '))^2=2(K_T(0)-K_T(\Theta )) =0,\end{aligned}$$

which implies that \(T(\theta ,\varphi ) = T(\theta ',\varphi ')\) with probability 1. Due to the isotropy of the random field \(T(\theta ,\varphi ),\) it is also true for any point on the sphere. \(\square \)

Theorem 2

The solution \(u(\eta , \theta , \varphi )\) of the initial random value problem (8), (18), (19) is given by the following random series

$$\begin{aligned} u(\eta ,\theta ,\varphi ) = \sum _{l=0}^{\infty } F_l(\eta ) \sum _{m=-l}^l a_{lm}Y_{lm}(\theta , \varphi ). \end{aligned}$$
(21)

The covariance function of \(u(\eta , \theta , \varphi )\) is given by

$$\begin{aligned} Cov(u(\eta , \theta , \varphi ), u(\eta ', \theta ', \varphi ')) =(4\pi )^{-1}\sum _{l=0}^\infty C_l (2l+1) F_l(\eta )F_l(\eta ')P_l(\cos \Theta ), \end{aligned}$$
(22)

where \(\eta ,\eta '\in [0,\eta _\infty ),\ \theta ,\theta '\in [0,\pi ], \ \varphi ,\varphi '\in [0,2\pi ),\) \(P_l(\cdot )\) is the lth Legendre polynomial and \(\Theta \) is the angular distance between the points \((\theta ,\varphi )\) and \((\theta ',\varphi ').\)

Remark 1

The convergence of the series in (21) is understood in the \(L_2(\Omega \times {\mathbb {S}}^2)\) sense. However, the series in (21) also converges almost surely, see [23, Section 2].

Proof

The solution of the initial value problem (8), (18), (19) is a spherical convolution of the function \({\widetilde{u}}(\cdot )\) obtained in Theorem 1 and the random field \(T(\theta ,\varphi ),\) provided that the corresponding Laplace series converges in the Hilbert space \(L_2(\Omega \times {\mathbb {S}}^2).\)

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)\) and have the Fourier–Laplace coefficients

$$\begin{aligned} a^{(i)}_{lm} = \int \limits _{{\mathbb {S}}^2}f_i(\theta ,\varphi )Y_{lm}^*(\theta , \varphi ) \sin \theta \textrm{d}\theta \textrm{d}\varphi , \quad i=1,2. \end{aligned}$$

The non-commutative spherical convolution of \(f_1(\cdot )\) and \(f_2(\cdot )\) is defined as the Laplace series (see [14])

$$\begin{aligned}{}[f_1*f_2](\theta , \varphi ) = \sum \limits _{l=0}^\infty \sum \limits _{m=-l}^l a_{lm}^*Y_{lm}(\theta , \varphi ) \end{aligned}$$
(23)

with the Fourier–Laplace coefficients given by

$$\begin{aligned} a_{lm} = \sqrt{\frac{4\pi }{2l+1}} a_{lm}^{(1)}a_{l0}^{(2)}, \end{aligned}$$

provided that the series in (23) converges in \(L_2(\Omega \times {\mathbb {S}}^2)\).

Thus, the random solution \(u(\eta , \theta , \varphi )\) of Eq. (8) with the initial conditions (18) and (19) can be written as a spherical random field with the following Laplace series representation

$$\begin{aligned} u(\eta , \theta , \varphi ) = [T*u](\theta ,\varphi ) = \sum _{l=0}^{\infty } F_l(\eta ) \sum _{m=-l}^l a_{lm}Y_{lm}(\theta , \varphi ). \end{aligned}$$

In the above series, the identity \(Y_{l0}^*({\textbf{0}}) = \sqrt{\frac{2\,l+1}{4\pi }}\) was applied to simplify the expression. By Lemma 3 and condition (20), the above series converges in \(L_2(\Omega \times {\mathbb {S}}^2).\)

By applying the addition formula (1) for the spherical harmonics, one obtains the covariance function of \(u(\eta ,\theta ,\varphi )\) in the form

$$\begin{aligned} Cov(u(\eta , \theta , \varphi ), u(\eta ', \theta ', \varphi ')) = \sum _{l=0}^\infty C_l F_l(\eta )F_l(\eta ')\sum _{m=-l}^lY_{lm}(\theta , \varphi ) Y^*_{lm}(\theta ', \varphi ')\\ = (4\pi )^{-1} \sum _{l=0}^\infty C_l (2l+1) F_l(\eta ) F_l(\eta ')P_l(\cos \Theta ).\end{aligned}$$

As \(|P_l(\cos (\Theta ))|\le 1,\) as \(l\rightarrow \infty ,\) it follows from (4) and Lemma 3 that the above series is convergent. \(\square \)

5 Properties of stochastic solution and its approximations

This section investigates time-smoothness properties of the solution and truncation errors of its approximation.

For \(L\in {\mathbb {N}},\) the following truncated series are used to approximate the solution of the random initial value problem in Theorem 2

$$\begin{aligned} u_L(\eta , \theta , \varphi ) = \sum _{l=0}^{L} F_l(\eta ) \sum _{m=-l}^l a_{lm}Y_{lm}(\theta , \varphi ). \end{aligned}$$
(24)

The next result provides the upper bounds for the corresponding approximation error.

Theorem 3

Let \(u(\eta , \theta , \varphi )\) be the solution of the initial value problem (8), (18), (19) and \(u_L(\eta , \theta , \varphi )\) be its approximation. Then, for \(\eta \in [0,\eta _\infty )\) the following bound for the truncation error holds true

$$\begin{aligned} \left| \left| u(\eta , \theta , \varphi ) - u_L(\eta , \theta , \varphi )\right| \right| _{L_2(\Omega \times {\mathbb {S}}^2)} \le C \left( \sum _{l=L+1}^{\infty } C_l (2l+1) \right) ^{1/2}, \end{aligned}$$

where the constant C does not depend on \(\eta .\)

Proof

The truncation error field \({\widehat{u}}_L(\eta , \theta , \varphi ) = u(\eta , \theta , \varphi ) - u_L(\eta , \theta , \varphi ), \ L \in {\mathbb {N}},\) is a centred Gaussian random field, i.e. \(E{\widehat{u}}_L(\eta , \theta , \varphi ) = 0\) for all \(L\in {\mathbb {N}},\ \theta \in [0,\pi ], \ \varphi \in [0,2\pi ),\) and \(\eta \in [0,\eta _\infty ).\) It follows from (21), (24) and the orthogonality of \(\{a_{lm}\}\) that

$$\begin{aligned} \left| \left| u(\eta , \theta , \varphi ) - u_L(\eta , \theta , \varphi )\right| \right| ^2_{L_2(\Omega \times S^2)} = \int \limits _{0}^{\pi }\int \limits _{0}^{2\pi } E\bigg ( \sum \limits _{l=L+1}^\infty F_l(\eta )\sum \limits _{m=-l}^l a_{lm} Y_{lm}(\theta , \varphi ) \bigg )^2 \sin \theta \textrm{d}\varphi \textrm{d}\theta .\end{aligned}$$

Then, due to the addition theorem for the spherical harmonics

$$\begin{aligned} ||u(\eta ,\theta ,\varphi )-u_L(\eta ,\theta ,\varphi )||^2_{L_2(\Omega \times {\mathbb {S}}^2)}=C\sum _{l=L+1}^{\infty } C_l (2l+1) F_l^2(\eta ).\end{aligned}$$

The statement of the theorem follows from Lemma 3. \(\square \)

Theorem 4

Let \(u(\eta , \theta , \varphi )\) be the solution of the initial value problem (8), (18), (19). If \(\sum \nolimits _{l=1}^\infty C_ll^3 <+\infty ,\) then for each \(\eta \in [0,\eta _\infty )\) and \(h\in (0,\eta _\infty -\eta )\)

$$\begin{aligned} ||u(\eta +h,\theta ,\varphi )-u(\eta ,\theta ,\varphi )||_{L_2({\mathbb {S}}^2\times \Omega )} \le Ch, \end{aligned}$$

where the constant C does not depend on \(\eta .\)

Proof

For \(h<\eta _\infty -\eta ,\) let us consider

$$\begin{aligned} ||u(\eta +h,\theta ,\varphi )-u(\eta ,\theta ,\varphi )||^2_{L_2(\Omega \times {\mathbb {S}}^2)} = \sum \limits _{l=1}^\infty C_l(2l+1) (F_l(\eta +h)-F_l(\eta ))^2.\end{aligned}$$

By (10) and applying the Cauchy inequality, one obtains that

$$\begin{aligned} ||u(\eta +h,\theta ,\varphi )-u(\eta ,\theta ,\varphi )||^2_{L_2(\Omega \times {\mathbb {S}}^2)} \le C \sum \limits _{l=1}^\infty C_l (2l+1) (A^2_l(\eta ,h)+B^2_l(\eta ,h)),\end{aligned}$$

where

$$\begin{aligned}A_l(\eta ,h):= K_1^{(l)}\big ( \eta _\infty - (\eta +h))^\nu J_\nu (z_l(\eta _\infty - (\eta +h))) - (\eta _\infty - \eta )^\nu J_\nu (z_l(\eta _\infty - \eta )) \big )\end{aligned}$$

and

$$\begin{aligned}B_l(\eta ,h):= K_2^{(l)}\big ( \eta _\infty - (\eta +h))^\nu Y_{\nu }(z_l(\eta _\infty - (\eta +h))) - (\eta _\infty - \eta )^\nu Y_{\nu }(z_l(\eta _\infty - \eta )) \big ).\end{aligned}$$

First, let us estimate \(|A_l(\eta ,h)|.\) Let \(f(x):=x^\nu J_{\nu }(x),\) then

$$\begin{aligned} A_l(\eta ,h) = K_1^{(l)}{z_l^{-\nu }}(f(z_l(\eta _\infty - (\eta +h)))-f(z_l(\eta _\infty - \eta ))).\end{aligned}$$

By the mean value theorem

$$\begin{aligned}|f(z_l(\eta _\infty - (\eta +h)))-f(z_l(\eta _\infty - \eta ))| \le z_l h \max _x|(x^\nu J_\nu (x))'|,\end{aligned}$$

where the maximum is taken over the interval \([z_l(\eta _\infty - (\eta +h)), z_l(\eta _\infty - \eta )].\) As \((x^\nu J_\nu (x))' = x^\nu J_{\nu -1}(x),\) one obtains

$$\begin{aligned}{} & {} {z_l^{-\nu }} |f(z_l(\eta _\infty - (\eta +h)))-f(z_l(\eta _\infty - \eta ))|\le {h}z_l^{-\nu +1} \\ {}{} & {} \quad \times \max _{\alpha \in [0,h]} |(z_l( \eta _\infty - (\eta +\alpha )))^\nu J_{\nu -1}(z_l(\eta _\infty - (\eta +\alpha )))|.\end{aligned}$$

As \(\sqrt{x}J_\nu (x)\) is a bounded function on \([0,\infty )\) and \(\nu >1,\) the next upper bound holds true

$$\begin{aligned} |A_l(\eta ,h)| \le |K_1^{(l)}| C {h}z_l^{-\nu +1} \max _{\alpha \in [0,h]} (z_l( \eta _\infty - (\eta +\alpha )))^{\nu -\frac{1}{2}} \le C|K_1^{(l)}| h \sqrt{z_l}.\end{aligned}$$

For \(B_l(\eta , h),\) by setting \(g(\eta ):= x^\nu Y_{\nu }(x)\) and using analogous calculations one obtains that

$$\begin{aligned} B_l(\eta ,h)\le & {} {z_l^{-\nu }}|K_2^{(l)}|\cdot |g(z_l(\eta _\infty - (\eta +h)))-g(z_l(\eta _\infty - \eta ))| \le {h}{z_l^{-\nu +1}}|K_2^{(l)}|\\ {}{} & {} \times \max _{\alpha \in [0,h]}| (z_l( \eta _\infty - (\eta +\alpha )))^\nu Y_{\nu -1}(z_l(\eta _\infty - (\eta +\alpha )))| \le {h}{z_l^{-\nu +1}}|K_2^{(l)}|\max _{x\in [0,z_l\eta _\infty ]}|x^\nu Y_{\nu -1}(x)|\\ {}\le & {} C{h}{z_l^{-\nu +1}}|K_2^{(l)}|(z_l\eta _\infty )^{\nu -\frac{1}{2}}\le C|K_2^{(l)}|hz_l^{\frac{1}{2}},\end{aligned}$$

as \(|x^\nu Y_{\nu -1}(x)|\) is bounded on \([0,A], \ A>0,\) and \(|x^\nu Y_{\nu -1}(x)|\le Cx^{\nu -\frac{1}{2}}\) when \(x\ge A.\)

Using the asymptotics of the functions \(J_\nu (\cdot )\) and \(Y_\nu (\cdot ),\) it is straightforward to verify that the constants \(K_1^{(l)}\) and \(K_1^{(2)}\) are bounded as \(|K_1^{(l)}|\le C\sqrt{z_l}\) and \(|K_2^{(l)}|\le C\sqrt{z_l}\) for all l. Thus, \(A^2_l(\eta ,h) \le Ch^2z^2_l\sim Ch^2l^2.\) Analogously, \(B^2_l(\eta ,h)\le Ch^2 z^2_l \sim Ch^2l^2\) which, by the assumption (4), implies the statement of the theorem. \(\square \)

6 On excursion probabilities

This section studies the excursion probabilities for the stochastic solution obtained in Sect. 4. The main approach used in this section is based on the metric entropy theory. This section also studies errors of approximations of extremes of \(u(\eta , \theta , \varphi )\) by extremes of the truncated field \(u_L(\eta , \theta , \varphi ).\) The approximation \(u_L(\eta , \theta , \varphi )\) converges to \(u(\eta , \theta , \varphi )\) in the space \(L_2(\Omega \times {\mathbb {S}}^2),\) as \(L\rightarrow \infty .\) However, in general, this type of convergence does not imply that the extremes of \(u(\eta , \theta , \varphi )\) can be effectively approximated using the extremes of \(u_L(\eta , \theta , \varphi ).\)

We first provide some results that will be used later to study properties of the distributions of extremes. In what follows T is a subset of \({\mathbb {R}}^n.\)

Theorem 5

([4, Theorem 2.1.1.]) Let \(\xi (x), \ x\in T,\) be a centred a.s. bounded Gaussian field, then

$$\begin{aligned} P\left( \sup \limits _{x\in T} \xi (x) - E\sup \limits _{x\in T} \xi (x) \ge y\right) \le \exp \left( {-\frac{y^2}{2\sigma _T^2}}\right) ,\quad y\ge 0, \end{aligned}$$

where \(\sigma _T^2:=\sup \limits _{x\in T}E\xi ^2(x),\) or, equivalently,

$$\begin{aligned} P\left( \sup \limits _{x\in T} \xi (x) \ge y \right) \le \exp \left( {-\frac{\left( y - E(\sup _{x\in T} \xi (x))\right) ^2}{2\sigma _T^2}}\right) , \quad y\ge E\sup \limits _{x\in T} \xi (x). \end{aligned}$$

Let us recall some results from the metric entropy theory. The canonical pseudometric generated by \(\xi (x), \ x\in T,\) is \( d(s,t) = \sqrt{ E(\xi (s) - \xi (t))^2 },\) \(s,t\in T.\)

Let \(B(t,\varepsilon ):=\{ s\in T: d(t,s)\le \varepsilon \}\) be a ball in T with the centre at \(t\in T\) and the radius \(\varepsilon .\) Then, the following Fernique’s inequality holds true [34, Theorem 2.5]

$$\begin{aligned} E\sup \limits _{x\in T} \xi (x) \le K\inf \limits _\mu \sup \limits _{x\in T} \int \limits _{0}^{\textrm{diam}(T)}\sqrt{\log \left( \frac{1}{\mu (B(x,\varepsilon ))}\right) }\textrm{d}\varepsilon ,\end{aligned}$$

where \(\mu \) is a probability measure on T,  \(\textrm{diam}(T)\) is the diameter of T in the metric \(d(\cdot )\) and K is a universal constant [35]. We will use this inequality to estimate \(E\sup \limits _{\theta ,\varphi }u(\eta , \theta , \varphi ).\) In the following results, \(\eta \in [0,\eta _\infty )\) is fixed.

By (22), the canonical pseudometric generated by \(u(\eta ,\theta ,\varphi )\) on the sphere is

$$\begin{aligned} d_\eta (\Theta ):= \left( E(u(\eta ,\theta ,\varphi ) - u(\eta ,\theta ',\varphi '))^2 \right) ^{\frac{1}{2}} = (\sqrt{2\pi })^{-1}\left( \sum _{l=0}^\infty C_l(2l+1) F_l^2(\eta )(1-P_l(\cos (\Theta ))) \right) ^\frac{1}{2}.\end{aligned}$$
(25)

Let us define “balls on the sphere” as \(B_\eta ((\theta ,\varphi ), \varepsilon ) = \{ (\theta ', \varphi '): d_\eta (\Theta )\le \varepsilon \},\) where \(\Theta \) is the angular distance between the points \((\theta ,\varphi )\) and \((\theta ',\varphi '),\) and also introduce the function

$$\begin{aligned} g_\eta (\varepsilon ) = \inf \{\Theta : \ d_\eta (\Theta ) \ge \varepsilon \}. \end{aligned}$$
(26)

For the simplicity of the exposition, let us denote the covariance function of \(u(\eta ,\theta ,\varphi )\) by

$$\begin{aligned} K_\eta (\Theta ):= Cov(u(\eta , \theta , \varphi ), u(\eta , \theta ', \varphi ')) = (4\pi )^{-1}\sum _{l=0}^\infty C_l (2l+1) F^2_l(\eta )P_l(\cos \Theta ).\end{aligned}$$

Note that by (25) \(K_\eta (\Theta )\) and \(d_\eta (\Theta )\) are linked by the relation \(K_\eta (0) - K_\eta (\Theta ) = d_\eta ^2(\Theta )/2. \)

Theorem 6

Let \(u(\eta , \theta , \varphi )\) be the solution of the initial value problem (8), (18), (19). If

$$\begin{aligned} \sum \limits _{l=0}^\infty C_ll^{2+\beta } <+ \infty \end{aligned}$$
(27)

with \(\beta >0,\) then

$$\begin{aligned} E\sup \limits _{\theta ,\varphi }u(\eta , \theta , \varphi ) \le K \int \limits _0^R \sqrt{\log \left( \frac{2}{1-\cos (g_\eta (\varepsilon ))}\right) }\textrm{d}\varepsilon , \end{aligned}$$
(28)

where \(R=\max _{\Theta \in [0,\pi ]} d_\eta (\Theta ),\) and the integral in (28) is finite.

Proof

By applying Fernique’s inequality, it follows that

$$\begin{aligned} E\sup \limits _{\theta ,\varphi }u(\eta , \theta , \varphi ) \le K\inf \limits _\mu \sup \limits _{\theta , \varphi } \int \limits _{0}^R\sqrt{\log \left( \frac{1}{\mu (B_\eta ((\theta , \varphi ),\varepsilon ))}\right) }\textrm{d}\varepsilon .\end{aligned}$$

The selection of \(\mu \) as the uniform distribution on the sphere results in the upper bound

$$\begin{aligned} E\sup \limits _{\theta ,\varphi }u(\eta , \theta , \varphi ) \le K\sup \limits _{\theta , \varphi } \int \limits _{0}^R\sqrt{\log \left( \frac{Leb({\mathbb {S}}^2)}{Leb(B_\eta ((\theta , \varphi ),\varepsilon ))}\right) }\textrm{d}\varepsilon .\end{aligned}$$

As \(u(\eta ,\theta ,\varphi )\) is an isotropic field, \(Leb(B_\eta ((\theta , \varphi ),\varepsilon ))\) is the same for all \(\theta , \varphi \) and can be replaced with \(Leb(B_\eta ({\textbf{0}},\varepsilon ))\). Note that \(B_\eta ({\textbf{0}},\varepsilon )\) is not necessarily a spherical cap as \(d_\eta (\Theta )\) is not necessarily strictly increasing in \(\Theta \). However, \(B_\eta ({\textbf{0}},\varepsilon )\) always contains a non-degenerate (not reduced to a single point) spherical cap. Indeed, \(u(\eta ,\theta ,\varphi )\) is \(L_2\)-continuous with respect to \(\theta \) and \(\varphi \) due to condition (27), see [23, Lemma 4.3]. Then, its decomposition (21) into spherical harmonics has the coefficients \(a_{lm}F_l(\eta ), \ l=0,1, \ldots ,\ m=-l, \ldots ,l,\) that are nonzero with probability 1, see Lemma 3. Hence, \(u(\eta ,\theta ,\varphi )\) is not a constant random variable over a set of \(\theta \) and \(\varphi ,\) and it follows from Lemma 5 that the spherical cap is non-degenerate.

The value of \( g_\eta (\varepsilon ) \in [0,\pi ]\) is the polar angle of the largest such cap. From the comparison of the corresponding areas, it follows that \(Leb(B_\eta ({\textbf{0}},\varepsilon )) \ge 2\pi (1-\cos (g_\eta (\varepsilon ))\) and

$$\begin{aligned} E\sup \limits _{\theta ,\varphi }u(\eta , \theta , \varphi ) \le K \int \limits _0^R \sqrt{\log \left( \frac{2}{1-\cos (g_\eta (\varepsilon ))}\right) }\textrm{d}\varepsilon . \end{aligned}$$

Let us provide conditions that guarantee the convergence of this integral for any fixed \(\eta \in [0,\eta _\infty ).\)

The Gaussian random field \(u(\eta ,\theta ,\varphi )\) is a.s. continuous and isotropic, and the corresponding spherical cap is non-degenerated. Due to Lemma 5, the continuous covariance function \(K_\eta (\Theta )\) is non-constant in any neighbourhood of the origin. Thus, by the definitions (25) and (26) \(g_\eta (\varepsilon ) \rightarrow 0,\) when \(\varepsilon \rightarrow 0.\)

As \(1-\cos (x) \sim x^2/2,\) \(x\rightarrow 0,\) it follows that

$$\begin{aligned}\log \left( \frac{2}{1-\cos (g_\eta (\varepsilon ))}\right) \sim -2\log (g_\eta (\varepsilon )), \ \ \ \varepsilon \rightarrow 0,\end{aligned}$$

which means that the integral in (28) is finite if and only if

$$\begin{aligned} \int \limits _0^R\sqrt{-\log (g_\eta (\varepsilon ))}\textrm{d}\varepsilon < \infty . \end{aligned}$$
(29)

By [23, Lemma 4.2], if \( \sum \limits _{l=0}^\infty C_ll^{2+\beta } < \infty \) for some \(0 < \beta \le 1,\) then

$$\begin{aligned} K_\eta (0) - K_\eta (\Theta ) = \frac{1}{2}d_\eta ^2(\Theta ) \le C \Theta ^{\beta +1}. \end{aligned}$$
(30)

As \(\log (1/\Theta ) \le 1/\Theta \) for sufficiently small \( \Theta ,\) it also holds true \(d_\eta (\Theta ) \le \frac{C}{(-\ln \Theta )^{(\beta +1)/2}}.\) Denoting \(\Theta = \exp ({-C\varepsilon ^{2\alpha }})\) for \(\alpha = -{1}/{(\beta +1)}\in (-1,-{1}/{2}],\) one gets \(d_\eta (\exp ({-C\varepsilon ^{2\alpha }}))\le \varepsilon .\) It follows from (26) that \(g_\eta (\varepsilon ) \ge \exp ({-C\varepsilon ^{2\alpha }})\) and \(-\log (g_\eta (\varepsilon )) \le C\varepsilon ^{2\alpha }.\) Thus, the integral (29) is finite, which finishes the proof of the theorem. \(\square \)

Now we are ready to state the main results of this section.

Theorem 7

For each fixed \(\eta \in [0,\eta _\infty )\), the following estimate holds true

$$\begin{aligned} P\left( \sup \limits _{\theta ,\varphi }u(\eta ,\theta ,\varphi ) > x \right) \le \exp \left( -\frac{\left( x-E(\sup _{\theta ,\varphi }u(\eta ,\theta ,\varphi ))\right) ^2}{2\sigma _\eta ^2}\right) ,\end{aligned}$$

where \(\sigma _\eta ^2 = (4\pi )^{-1}\sum _{l=0}^\infty C_l (2\,l+1) F^2_l(\eta ).\)

If \( \sum \nolimits _{l=0}^\infty C_ll^{2+\beta } <+ \infty \) with \(\beta > 0,\) then

$$\begin{aligned} P\left( \sup \limits _{\theta ,\varphi }u(\eta ,\theta ,\varphi ) > x \right) \le \exp \left( -\frac{\left( x- K_1\right) ^2}{2\sigma _{\eta ,L}^2} \right) \end{aligned}$$

for \(x> K_1:=K \int \limits _0^R \sqrt{\log \left( \frac{2}{1-\cos (g_{\eta }(\varepsilon ))}\right) }\textrm{d}\varepsilon ,\) where the latter integral is finite.

Proof

The statement is a direct corollary of Theorems 5 and 6 as the application of the estimate for \(E\sup \limits _{\theta ,\varphi }u(\eta ,\theta ,\varphi )\) given by (28) can only increase the upper bound when \(x > K_1.\) \(\square \)

Now, let us estimate of probabilities of large deviations between extremes of the stochastic solution \(u(\eta , \theta , \varphi )\) and its approximation \(u_L(\eta , \theta , \varphi ).\) The truncation error field \({\widehat{u}}_L(\eta , \theta , \varphi ) = {u}(\eta , \theta , \varphi ) - u_L(\eta , \theta , \varphi ), \ L \in {\mathbb {N}},\) is a centred Gaussian random field. Similar to the considered results, it generates the canonic pseudometric \(d_{L,\eta }(\Theta )\) on \({\mathbb {S}}^2\) and the function \(g_{L,\eta }(\varepsilon )\) is associated with \(d_{L,\eta }(\cdot ),\) which is given by \(g_{L,\eta }(\varepsilon ):=\min \{\Theta : d_{L,\eta }(\Theta ) \ge \varepsilon \}.\) Analogously to the previous results, one obtains:

Corollary 1

For each fixed \(\eta \in [0,\eta _\infty )\) the following estimate holds true

$$\begin{aligned} P\left( \sup _{\theta ,\varphi }|{\widehat{u}}_L(\eta ,\theta ,\varphi )| > x \right) \le 2\exp \left( -\frac{\left( x-E\left( \sup _{\theta ,\varphi }{\widehat{u}}_L(\eta ,\theta ,\varphi )\right) \right) ^2}{2\sigma _{\eta ,L}^2}\right) \end{aligned}$$
(31)

where \(\sigma _{\eta ,L}^2 = (4\pi )^{-1}\sum _{l=L+1}^\infty C_l (2\,l+1) F^2_l(\eta ).\)

If there exists \(\beta > 0,\) such that

$$\begin{aligned} \sum \limits _{l=0}^\infty C_ll^{2+\beta } <+ \infty ,\end{aligned}$$
(32)

then

$$\begin{aligned}P\left( \sup \limits _{\theta ,\varphi }|{\widehat{u}}_L(\eta ,\theta ,\varphi )| > x \right) \le 2\exp \left( -\frac{\left( x- K_2\right) ^2}{2\sigma _{\eta ,L}^2} \right) \end{aligned}$$

for \(x> K_2:=K \int \limits _0^R \sqrt{\log \left( \frac{2}{1-\cos (g_{L,\eta }(\varepsilon ))}\right) }\textrm{d}\varepsilon ,\) where the latter integral is finite.

Proof

As \({\widehat{u}}_L(\eta ,\theta , \varphi )\) is the centred Gaussian random field, it holds for \(x>0\) that

$$\begin{aligned} P\left( \sup _{\theta ,\varphi }|{\widehat{u}}_L(\eta ,\theta , \varphi )|> x \right) \le 2P\left( \sup _{\theta ,\varphi }{\widehat{u}}_L(\eta ,\theta , \varphi ) > x \right) . \end{aligned}$$

As

$$\begin{aligned} \sigma _{\eta ,L}^2 = E {\widehat{u}}^2_L(\eta ,\theta , \varphi ) = (4\pi )^{-1}\sum _{l=L+1}^\infty C_l (2l+1) F^2_l(\eta ),\end{aligned}$$

the inequality (31) is a direct corollary of Theorem 5.

The random field \({\widehat{u}}_L(\eta ,\theta , \varphi )\) is \(L_2\)-continuous in \(\theta \) and \(\varphi \) by condition (32). By applying Fernique’s inequality with the uniform distribution \(\mu \), one obtains

$$\begin{aligned} E\sup \limits _{\theta ,\varphi }{\widehat{u}}_L(\eta , \theta , \varphi ) \le K \int \limits _0^R \sqrt{\log \left( \frac{2}{1-\cos (g_{L,\eta }(\varepsilon ))}\right) }d\varepsilon . \end{aligned}$$

The convergence of the above integral can be shown by using steps analogous to the proof of Theorem 6. The application of the above estimate to (31) completes the proof. \(\square \)

7 Numerical studies

This section presents numerical studies of the solution \(u(\eta ,\theta ,\varphi )\) of the initial value problem (8), (18) and (19) and its truncated approximation \(u_L(\eta , \theta , \varphi )\) from Sect. 5.

Numerical computations were performed using the software R version 4.3.1 and Python version 3.11.5. The HEALPix representation (http://healpix.sourceforge.net) and the Python package “healpy” were used for computations. The Python package “pyshtools” was used to simulate spherical random fields. The R package “rcosmo”, see [16, 17], was used for visualizations. The R and Python code are freely available in the folder “Research materials” from the website https://sites.google.com/site/olenkoandriy/. For the numerical studies, we use a map of CMB intensities obtained by the mission Planck. These data are available in [20] as well as its angular power spectrum in [21].

7.1 Functions \(F_l(\eta )\)

The evolution in time of the solution \(u(\eta ,\theta ,\varphi )\) is governed by the functions \(F_l(\eta ),\) \(\eta \in [0,\eta _\infty ),\) \(l=0,1 \ldots ,\) given by (10). These functions are the multiplication factors defining the change of the angular power spectrum of the solution with time \(\eta \). Lemma 2 shows that for a fixed \(\eta \) the multiplication factors have a wave structure for large l. Figure 1a depicts multiplication factors \(F_\eta (l)\) as functions of the argument l for the fixed time moments \(\eta =0.001\) and \(\eta =0.002.\) One can see that \(F_\eta (l),\) as the functions of the argument l,  form waves whose periods decrease in time. On the other hand, Fig. 1b displays \(F_\eta (l)\) as functions of the argument \(\eta \) for the fixed values of indices \(l=3\) and \(l=10.\) In this case, the functions \(F_\eta (l)\) form waves whose amplitudes decrease as they propagate. Figure 1b also demonstrates that their periods decrease with l.

Fig. 1
figure 1

\(F_\eta (l)\) as functions of arguments l and \(\eta \)

Full size image

7.2 Evolution of the solution

This subsection provides numerical simulations for the evolution in space and time according to the model (8), (18) and (19), and studies spatio-temporal dependencies and evolution of the angular spectrum. For illustrative purposes, the constants \(c,D,r,\eta _\infty \) were set equal to 1.

The initial condition field consisting of the measurements of CMB intensities was used. The structure of the CMB data is quite complex with various hot and cold areas located close to each other. Figure 2a shows the map of CMB. From the mathematical point of view, these values of the CMB intensities form a single realization of a spherical isotropic Gaussian random field. Thus, its spherical harmonics representation (2) can be derived. The Python’s package “healpy” was used to compute the coefficients \(a_{lm}\) from the map of CMB intensities. In the following numerical examples the values \(a_{lm}\) for \(l=0,1, \ldots ,2500,\) and \(m=-l, \ldots ,l\) and the corresponding initial condition \(u(0,\theta ,\varphi ) =\sum _{l=0}^{2500}\sum _{m=-l}^l a_{lm} Y_{lm}(\theta , \varphi )\) were used.

The representation (21) was used to obtain the solution \(u(\eta , \theta , \varphi )\) for \(\eta =0.001.\) From the visualization of \(u(0.001,\theta ,\varphi )\) in Fig. 2b, one can see that the solution becomes smoother and large deviations (both positive and negative) from the overall mean are decreasing in time.

Fig. 2
figure 2

Realizations of \(u(\eta ,\theta ,\varphi )\)

Full size image
Fig. 3
figure 3

Spatial dependencies and coefficients \(a_{lm}\) for \(\eta =0\) and \(\eta =0.001\)

Full size image

Figure 3a demonstrates how spatial dependencies change with increasing angular distance between spherical locations. The spatial correlation function decreases rapidly with the increase in the angular distance \(\Theta .\) Figure 3a, b shows the change of spatial correlations for \(\eta =0\) and \(\eta =0.001.\) Figure 3c, d plot values of the spherical harmonic coefficients \(a_{lm}\) of the realizations in Fig. 2 for levels up to 200.

Fig. 4
figure 4

Angular spectra of \(u(\eta ,\theta ,\varphi )\) for \(\eta =0\) and \(\eta =0.001\)

Full size image

We also employed the angular power spectrum \(C_l\) of CMB provided in [21]. It is depicted in Fig. 4 when \(\eta =0\). Then, the spherical angular power spectrum for the case of \(\eta =0.001\) was computed. The plots of spherical angular power spectra for \(\eta =0\) and \(\eta =0.001\) show that they remain almost identical for small values of l and flatten out when l and \(\eta \) increase. This observation aligns with the obtained theoretical findings and cosmological theories positing that higher multipoles undergo more rapid changes.

To illustrate the structure of space–time dependencies, we also produced a 3D plot showing normalized covariances (divided by \(Eu^2(0,0,0))\) as a function of angular distance \(\Theta \) and time \(\eta .\) In Fig. 5, when \(\Theta \) and \(\eta \) increase, the normalized covariance function rapidly decreases for small values of \(\Theta \) and \(\eta ;\) whereas, for other values, the decay of covariances is rather slow. One can see that the variance of the solution \(Eu^2(\eta ,0,0)\) is decaying with \(\eta \) rapidly for small values of \(\eta .\)

Fig. 5
figure 5

Spatio-temporal covariances of \(u(\eta ,\theta ,\varphi ).\)

Full size image

7.3 Excursion probabilities

To analyse the excursion probabilities, we used the solutions with the initial value field possessing the power spectrum shown in Fig. 4 for the case \(\eta =0.\) By setting the value \(\eta =0.001\) and by using the relationships (25) and (26), one can compute the canonical pseudometric \(d_{0.001}(\Theta )\) generated by \(u(0.001,\theta ,\varphi )\) on the sphere, along with the associated function \(g_{0.001}(\varepsilon ).\) The results of these computations are illustrated in Fig. 6. The plots demonstrate a rapid increase of the function \(d_{0.001}(\Theta )\) in the neighbourhood of the origin. This observed behaviour of the canonical pseudometric can be attributed to the rapid change of the spatial covariance function \(K_{0.001}(\Theta )\) at the origin, see Fig. 3a and consult (25) on their relations.

Fig. 6
figure 6

Canonical pseudometric \(d_\eta (\Theta )\) and the function \(g_{\eta }(\varepsilon )\) for \(\eta =0.001\)

Full size image

By using the Python package “pyshtools” and the values \(C_l\) of the CMB angular power spectrum, we simulate 300 realizations of spherical Gaussian random fields. Then, those realizations were used as the initial conditions for the problem (8), (18) and (19). By decomposing each of these 300 realizations into a series of spherical harmonics and by applying the formula (21) with \(\eta =0.001\), we obtained 300 realizations of the solution for \(\eta =0.001.\)

Fig. 7
figure 7

Sample distribution of \(\sup _{\theta ,\varphi }u(0.001, \theta ,\varphi )\)

Full size image
Fig. 8
figure 8

Sample probabilities \(P(\sup _{\theta ,\varphi }u(0.001, \theta ,\varphi )>x)\) and theoretical bounds

Full size image

Figure 7 shows the sample distribution of \(\sup _{\theta ,\varphi }u(0.001, \theta ,\varphi )\) obtained by using the realizations of the solution. It is known that the distribution of the supremum of a Gaussian process on an interval is positively skewed as its median is less than or equal to its mean. The similar behaviour for the case of spherical Gaussian random fields is observed in Fig. 7. Figure 8 provides the estimated excursion probabilities \(P(\sup _{\theta ,\varphi }u(0.001, \theta ,\varphi )>x)\) and the corresponding upper bound given in Theorem 7. The observed pattern is analogous for this type of estimators. For small values of x,  the upper bound appears conservative; however, it rapidly approaches the excursion probability as x increases. Figure 8 illustrates the areas for potential improving of this bound.

8 Conclusions and future work

The paper investigated spherical diffusion within an expanding space–time framework. It derived solutions to the deterministic and stochastic diffusion equations, specifically focusing on the exponential growth case. The findings offer insights into the asymptotic and local characteristics of the solutions. Additionally, the paper explored the extremal properties of these solutions. Given the validity of Borell–TIS and Fernique’s inequalities for general sets, a similar methodology can be employed to investigate extremes of solutions on more complex manifolds. Applications of the theoretical findings were demonstrated through simulation studies by modelling evolutionary scenarios and examining properties of the obtained estimates with the CMB intensities as the initial condition. In future research, it would be interesting to further study properties of this model and the corresponding solutions and to extend the developed approach to other SPDEs.