On approximation for time-fractional stochastic diffusion equations on the unit sphere

This paper develops a two-stage stochastic model to investigate evolution of random fields on the unit sphere $\bS^2$ in $\R^3$. The model is defined by a time-fractional stochastic diffusion equation on $\bS^2$ governed by a diffusion operator with the time-fractional derivative defined in the Riemann-Liouville sense. In the first stage, the model is characterized by a homogeneous problem with an isotropic Gaussian random field on $\bS^2$ as an initial condition. In the second stage, the model becomes an inhomogeneous problem driven by a time-delayed Brownian motion on $\bS^2$. The solution to the model is given in the form of an expansion in terms of complex spherical harmonics. An approximation to the solution is given by truncating the expansion of the solution at degree $L\geq1$. The rate of convergence of the truncation errors as a function of $L$ and the mean square errors as a function of time are also derived. It is shown that the convergence rates depend not only on the decay of the angular power spectrum of the driving noise and the initial condition, but also on the order of the fractional derivative. We study sample properties of the stochastic solution and show that the solution is an isotropic H\"{o}lder continuous random field. Numerical examples and simulations inspired by the cosmic microwave background (CMB) are given to illustrate the theoretical findings.


Introduction
Let  2 ⊂ ℝ 3 be the 2-dimensional unit sphere in the Euclidean space ℝ 3 .In this paper, we develop a two-stage stochastic model defined by time-fractional stochastic diffusion equations on  2 .The stochastic model can be used to model the evolution of spatio-temporal stochastic systems such as climate changes and the density fluctuations in the primeval universe (see [2,4,5,28]).For general stochastic diffusion models (or Itô diffusion) we refer the reader to [27]).
The model is constructed as follows.Let  be a strongly isotropic Gaussian random field on  2 and let  (),  ∈ (0, ∞), be a random field on  2 , which solves the following time-fractional diffusion equation: with the initial condition  (0) = , where Δ  2 is the Laplace-Beltrami operator,  is a constant satisfying 0 <  ≤ 1,    is the Riemann-Liouville time fractional derivative of order 0 ≤  < 1 defined, see [31], as and   ,  > 0, is a time-delayed Brownian motion defined in  2 (Ω ×  2 ) (the product space of the probability space Ω and the sphere  2 ).The structure of the noise   will be given in Section 2.3.
In words, the model in (1) describes an initial isotropic random field evolving by a time-fractional diffusion process, with no external input until time , at which time external noise is switched on.
Fractional stochastic diffusion equations on the sphere Note that in the model (1), it is assumed that the random field  is independent of the noise   .Since equation ( 1) is linear, it is easy to see that  can be decomposed as  () ∶=   () +   (),  ∈ (0, ∞), where   satisfies the homogeneous equation with the initial condition   (0) = , and   satisfies the inhomogeneous equation with the condition   () = 0,  ∈ (0, ], from which it follows that   () is zero until time .
We write the equation ( 5) in the integral form, for  > , as Note that in (1) there are two sources of randomness, namely the initial condition  and the driving noise   .The split of the stochastic solution  in equation ( 3) into two independent parts (due to linearity) separates the two sources of randomness: the homogeneous solution   is random only through the initial condition, while the inhomogeneous solution   is random only through the time-delayed Brownian motion   .We shall derive the stochastic solution  to the equation (1), which represents a random field defined in  2 (Ω× 2 ).The solution  is given in the form of an expansion in terms of complex orthonormal spherical harmonics { , ∶  ∈ ℕ 0 ,  = −, … , }, i.e.,  (, ) ∶= ∞ ∑ =0  ∑ =− Û, (, ) , ,  ∈ (0, ∞), (7) where the random variables Û, represent the Fourier coefficients of  and are given by Û, (, ) ∶= ∫  2  (, , ) , ()(),  ∈ (0, ∞), where  , denotes the complex conjugate of  , and  is the normalized Riemann surface measure on the sphere  2 (details are given in Section 2).Given that the initial field  is a strongly isotropic Gaussian random field on  2 , we shall show that the stochastic solution  is also a strongly isotropic Gaussian random field on  2 .
We truncate the expansion (7)  We investigate the approximate solution   and its rate of convergence as  → ∞.We shall show that the convergence rate (in  2 -norm on  2 (Ω ×  2 )) depends not only on the decay of the power spectrum of the driving noise   and the initial condition , but also on the order of the fractional derivative .Also, we show that the stochastic solution  (, ),  > , evolves continuously with time by deriving an upper bound in the  2 -norm for its temporal increments from time  to  + ℎ of the form ()ℎ 1 2 , ℎ > 0,  ∈ (, ∞), where the continuous function (⋅) does not depend on ℎ.Finally, under some conditions, we show the existence of a locally Hölder continuous modification of  .
Note that the proposed model, defined in (1), describes the evolution (in time) of a two-stage stochastic system.The homogeneous equation given in (4) determines the evolution for the system with no external input while the inhomogeneous equation (5) provides the perturbation produced by the Brownian motion starting at time .The cosmic microwave background (CMB) provides motivation for considering multi-stage stochastic systems.Roughly speaking, CMB is electromagnetic radiation that has been emitted across the cosmos since ionised atoms and electrons recombined roughly 370,000 years after the Big Bang [9,29,30].Since CMB has passed through several stages called formation epochs or phases (such as Planck, Grand Unification, Inflation and recombination), it can be considered as an example of a multi-stage stochastic system.
The literature shows a variety of studies where models using stochastic partial differential equations were developed, see for example [4,6,18,21,22,25].However, little attention has been paid to two-stage stochastic models of the kind considered here.A fractional SPDE, governed by a fractional derivative in time and a fractional diffusion operator in space, was developed in [1] to model evolution of tangent vector random fields on the unit sphere.In a recent paper, Anh et.al.[2] considered a two-stage stochastic model defined by SPDEs on the sphere governed by a fractional diffusion operator in space and a fractional Brownian motion as driving noise.The case where a random initial condition given by a fractional stochastic Cauchy problem was considered.They showed that the truncation errors of the stochastic solution have, in the  2 −norm, the convergence rates    − ,  > 1, under the assumption that the variances of the driving noise satisfy some smoothness condition.Also, as a function of time , it was demonstrated that the "constant"   blows up when  is sufficiently small.In a recent article [16], a time-fractional stochastic heat equation driven by time-space white noise in ℝ  was considered.The model is defined in the distribution sense, and an explicit solution is derived within the space of tempered distributions.The stochastic solution is demonstrated to exhibit mild behavior exclusively when the dimension is either  = 1 or  = 2, and the order of the fractional derivative falls within the range  ∈ (1, 2), while for  < 1, the solution is not mild for any .
The paper is structured as follows.Section 2 presents necessary material from the theory of functions on the sphere  2 , spherical harmonics, Gaussian random fields,  2 ( 2 )-valued time-delayed Brownian motions, and some necessary tools from the theory of stochastic integrals.Section 3 derives properties of the solution to the homogeneous equation (4).In Section 4 we study the solution to the inhomogeneous equation (5).In Section 5 we give results for the combined equation (1).Section 6 derives an approximation to the solution of equation (1), the rate of convergence for the truncation errors.In Section 7 we study the temporal increments, in the  2 -norm, of the stochastic solution of (1) and its sample Hölder continuity property.In Section 8 we provide some numerical examples to explain the theoretical findings.In particular, Section 8.1 illustrates the evolution of the stochastic solution of (1) using simulated data inspired by the CMB map.Section 8.2 explains the convergence rates of the truncation errors, in the  2 -norm, of the stochastic solution of (1).Finally, Section 8.3 explores the convergence rates of the temporal increments, in the  2 -norm, of the solution of (1).
The basis  , and the Legendre polynomial   satisfy the addition theorem (see [26]), that is for ,  ∈  2 , there holds The Laplace-Beltrami operator Δ  2 (or the spherical Laplacian) on the sphere  2 at  is given in terms of the spherical coordinates by (see [8,26]) It is well-known that the spherical harmonic functions { , ∶  ∈ ℕ 0 ,  = −, … , } are the eigenfunctions of the negative Laplace-Beltrami operator −Δ  2 on the sphere  2 with eigenvalues that is for An arbitrary complex-valued function  ∈  2 ( 2 ) can be expanded in terms of a Fourier-Laplace series, in the  2 ( 2 ) sense, The f, , for  ∈ ℕ 0 ,  = −, … , , are known as the Fourier coefficients for the function  under the Fourier basis  , .By Parseval's theorem, for  ∈  2 ( 2 ), there holds Note that by the properties of the spherical harmonic functions, the Fourier-Laplace series (11) of any real-valued  ∈  2 ( 2 ) takes the form )) .

Isotropic random fields on the sphere
This subsection introduces isotropic Gaussian random fields on the sphere  2 and their expansions in terms of complex spherical harmonics [12,21].
In this paper, we consider a particular class of random fields called isotropic random fields.Following [24], the random field  is strongly isotropic if for any  ∈ ℕ and for all sets of  points  1 , … ,   ∈  2 , and for any rotation  ∈ (3), the joint distributions of ( 1 ), … , (  ) and ( 1 ), … , (  ) coincide.
A random field  on  2 is Gaussian if for each  ∈ ℕ = {1, 2, 3, … } and each  points  1 , … ,   ∈  2 , the vector (( 1 ), … , (  )) has a multivariate Gaussian distribution.Note that a Gaussian random field  is strongly isotropic if and only if it is 2-weakly isotropic (see [24]).In this paper, we consider random fields that are Gaussian and strongly isotropic.
• The elements of  ′ with  > 0 have the property that Re ξ, and Im ξ, are independent and  (0,   ∕2) distributed.
Note that it is easy to show that   (),  ≥  ≥ 0, is a centered Gaussian process on [, ∞) satisfying Putting  =  in (17) we have The case when  = 0,  0 (),  ≥ 0, represents the standard Brownian motion.We will for brevity write  0 () as () if no confusion arises.We define a time-delayed Brownian motion on the unit sphere PROOF.By definition, In what follows we will assume that the condition (19) remains valid.
Using the above setting, we have the following proposition.Its proof is similar to the one from [7,Proposition 4.3].
Note that the sequence {  ∶  ∈ ℕ 0 } is called the angular power spectrum of   .
In subsequent sections, some results are required for establishing the pathwise solution of (5).In particular, a formal justification for taking the fractional derivative of stochastic integrals defined by ( 29) is required.Since there is currently no existing result in the literature on this matter, we embark on proving the following theorem.
and by Proposition 2.5, the stochastic integral is well-defined.Let us show that the stochastic integral is well-defined.Using Itô's isometry and Lemma 3.10 in [31] we obtain where the last second step used the relations (4.3.1) and (5.1.14)in [13] and the upper bound (26).By the relation (39) we obtain Using equations (2.88) and (1.100) in [31], Thus, equation (41) becomes which completes the proof.

Solution of the homogeneous problem
Consider a random field   () ∈  2 (Ω ×  2 ),  ∈ (0, ∞).We say that the field   satisfies the equation ( 4), in the  2 (Ω ×  2 ) sense, if for a given  > 0, there holds This section derives the solution   (),  ∈ (0, ∞), to the homogeneous equation ( 4) under the random initial condition   (0) = , where  is a centered, strongly isotropic Gaussian random field on  2 .Note that since   ∈  2 (Ω ×  2 ), it follows that the field   takes the series expansion where the random coefficients Û  , ,  ∈ ℕ 0 ,  = −, … , , are given by By multiplying both sides of (4) by  , and integrating over  2 , we obtain, with the help of (10), the following set of ordinary differential equations We denote by Ũ  , () the Laplace transform of Û  , with respect to , i.e., assuming that the integral converges.The condition for its convergence will be given in Proposition 3.1.Since 0), then by taking the Laplace transform of both sides of (42), we get Now solving for Ũ  , (), we arrive at Recall that (see [31], equation 1.80) where   (⋅) is the Mittag-Leffler function (see equation ( 25)).By taking the inverse Laplace transform in (44), the solution of (4) becomes In order to prove that   is the solution of (4) we need to prove a truncation field    of   is a solution and then pass it to the limit as the truncation degree  → ∞.Also, we need the uniform convergence of      to     in some sense.We begin with the following lemma.
Let    () be defined as where ξ, are the Fourier coefficients of .
Let   be defined as Using (54) with the Parseval formula, (13) and the upper bound (33) then for  <  we have For a given  > 0, using the condition (15), there is  0 independent of  such that for ,  ≥  0 , the RHS of (55) can be smaller than .We then let which is convergent in  2 (Ω ×  2 ).
Proposition 3.1.Let the angular power spectrum {  ∶  ∈ ℕ 0 } of the isotropic Gaussian random field  satisfy assumption (15).Then the random field   defined by (45) satisfies (4), in the  2 (Ω ×  2 ) sense, under the initial condition   (0) = .In particular, for a given  0 > 0, we have PROOF.Let us fix  ≥ 1, then we have Using Lemma 3.10 in [13] we have Combining the above two equations, then for any fixed  ∈ Ω we obtain It follows from (57) that By the triangle inequality ] .
By Lemma 3.2, we can choose  so that By Lemma 3.4, we can choose  so that So the proposition is proved.
The following result shows that the solution   to the equation ( 4) is a centered, 2-weakly isotropic Gaussian random field on  2 .Here we assume that the angular power spectrum {  ∶  ∈ ℕ 0 } of  decays algebraically with order  1 > 2, i.e., there exist constants D, C > 0 such that (60) For  1 > 2 and C > 0, we let Proposition 3.2.Let the field   , given in (45), be the solution to the equation (4) under the initial condition   (0) = , where  is a centered, 2-weakly isotropic Gaussian random field on  2 .Let {  ∶  ∈ ℕ 0 }, the angular power spectrum of , satisfy (60).For a fixed  ∈ (0, ∞),   () is a centred, 2-weakly isotropic Gaussian random field on  2 , and its random coefficients where   is given in (9) and   ′ is the Kronecker delta function.
For  ∈ ℕ, let where the third step used Lemma 3.10 in [13].
PROOF.For  ≥ ,  < , by Parseval's formula and Itô's isometry, we have 2  by using ( 13) where the third step used (26).For a given  > 0, using the condition (19), there is  0 independent of  such that for ,  ≥  0 , the RHS of (78) is smaller than .So {   } is a Cauchy sequence in  2 (Ω ×  2 ) and hence it is convergent.We define the limit of    as  → ∞ as in (77).Note that by the condition ( 19) we obtain which guarantees that   is well-defined in the  2 (Ω ×  2 ) sense.

Solution of the time-fractional diffusion equation
This section demonstrates the solution  (),  ∈ (0, ∞), to the time-fractional diffusion equation given in (1).

Approximation to the solution
The results in the previous sections give a series representation of the solution  (),  ∈ (0, ∞), of the timefractional diffusion equation (1).This section provides an approximation to the solution  () by truncating its expansion at a degree (truncation level)  ∈ ℕ.Then we give an upper bound for the approximation error of   .Definition 6.1.The approximation   of degree  ∈ ℕ to the solution  given in (92) is defined as where    () is given by (94) and for  ≤ ,    () = 0, while for  > ,    () is given by (95).Proposition 6.1.Let {  ∶  ∈ ℕ 0 }, the angular power spectrum of , satisfy (60).Let    (),  ∈ (0, ∞), be defined as , where    is given in (94).Then, the following estimates hold true: where C 1 is given by (61).
b) for  > where    (⋅) is defined as PROOF.By (94) we can write ] .
By the properties of the spherical harmonics we have and hence 99) can be bounded by Now, let us consider part b.Since   (−    ) is positive and decreasing for  > 0, then for  >  − 1   , using ( 33), (99) can be bounded by thus completing the proof.
In the next proposition, we derive an upper bound for the approximation errors of the truncated solution    (),  > .
To prove ii, note that by using the upper bound (31) with the estimates (34), (35), and (36), we obtain where   is given in (27).By using (105) with the estimate (106) and the properties of   given in (84), we get where   ( 2 ) is given in (86) and    (⋅) is given in (102), thus completing the proof.
Remark 6.1.In view of Theorem 6.1 and the functions    () and    () given by (97) and (102) respectively, it is easy to see that for  ∈ ( 12 , 1], both functions    () and    () are bounded when  → ∞ and hence the truncation errors   () are bounded when  → ∞ (see case III).Also, it is worth noting that when  ∈ ( 12 , 1], the truncation errors   () as functions of  are bounded as well when  → 0 (see case I), which improves the upper bound for the truncation errors obtained in [2] when  = 1 and  is sufficiently small.

Temporal increments and sample Hölder continuity of solution
In this section we study some properties of the stochastic solution of equation (1).In particular, we derive an upper bound, in  2 -norm, of the temporal increments of the stochastic solution  (),  ∈ (, ∞), to equation (1).Then we study the sample properties of the stochastic solution  .Under some conditions, we show the existence of a locally Hölder continuous modification of the random field  .We demonstrate how the ℙ-almost sure Hölder continuity of the solution  depends on the decay of the angular power spectra of the initial condition  and the driving noise   .Theorem 7.1.Let  be the solution (92) to the equation (1).Let  ℎ , ℎ > 0, be defined as ,  > . (109) Then, for ℎ > 0, there holds where (),  > , is given by ,  > 0, C 1 , and Ã 2 are defined in (61) and (85) respectively.

Numerical studies
In this section, we present some numerical examples for the solution  of equation (1).In particular, we show the evolution of the stochastic solution  of equation (1) using simulated data inspired by the CMB map.Also, we explain the convergence rates of the truncation errors of the approximations   (),  ∈ (0, ∞), and the  2 -norm of the temporal increments of the stochastic solution  (),  ∈ (, ∞), to the equation (1).

Evolution of the solution
In this subsection, we illustrate the evolution of the stochastic solution  of the equation (1) using simulated data inspired by the CMB map.
First, in Figure 1 we display a realization of the truncated initial condition at time  = 0 of degree  = 600 (i.e. =  600 (0)) whose angular power spectrum   has the form given by (60) with D = C = 1,  1 = 2.3.The Python HEALPy package (see [14]) was used to generate such a realization.Then, we use the obtained initial realization  600 (0) to generate two realizations for truncated solution  600 () given by (96) at different times, namely  =  and  = 10 with  = 10 −5 for the case  = 0.5,  1 = 2.3, and  2 = 2.5.We use the equation (84) with K = Ã = 10 4 to compute the angular power spectrum   for the time-delayed Brownian motion   .The Fourier coefficients Û, () for the solution   () are computed, using (93) with ( 68) and ( 90 where  , ∶=  (1)  , − i (2)  , and  ,, () ∶=  (1)  ,, () − i (2)  ,, (),  > 0. Note that for each realization and given time  > , we simulate the stochastic integrals  () ,, ( − ),  = 1, 2, using Proposition 2.5, as independent, normally distributed random variables with mean zero and variances  2 ,−, where  2 ,, is given by (28).Figures 2a and 2b, respectively, show two realizations of the truncated solutions   () and    () with  = 600, using the initial realization obtained in Figure 1, at times  =  and  = 10 with  = 10 −5 .Figure 3a shows a realization of the inhomogeneous solution    () at time  = 10 with  = 10 −5 and  = 0.5, while Figure 3b shows a realization of the combined solution   () using the realizations obtained in Figures 2b and 3a   The Figures 2a and 2b reflect the time evolution of the truncated homogeneous stochastic solution   600 () at time  =  and  = 10 respectively.In this case, the randomness comes only from the initial condition.Figure 3a reflects the time evolution of the truncated inhomogeneous stochastic solution   600 () at time  = 10 (in this case the randomness comes only from   ). Figure 3b reflects the evolution of the same realization obtained in Figure 2b at  = 10 but for which the noise   is switched on (in this case the randomness comes from both the initial condition and   , i.e.  600 (10) =   600 (10) +   600 (10)).

Simulations of truncation errors of solution
This subsection presents some numerical examples for the truncation errors of the combined stochastic solution  of the equation (1).In particular, the numerical examples explain the convergence rates of the truncation errors   (), given by (108), of the stochastic solution of (1).To produce numerical results, we use the angular power spectra   and   given by ( 60) and ( 84

Simulations of temporal increments
This subsection presents a numerical example for the stochastic solution  of the equation (1) in time.In particular, the numerical example explains the  2 −norm temporal increments  ℎ (), given by (109), of equation (1).
Figure 4  shows numerical errors  , L() of N = 100 realizations and the corresponding theoretical errors for different values of  and  with  1 = 2.3 and  2 = 2.5 for the case  = 0.5.In particular, in Figure4awe illustrate the case I (0 <  ≤  −2  ) in Theorem 6.1 with  = 10 −12 .In Figure4bwe illustrate the case II with  =  + −2   for  = 10 −5 , while Figure4cillustrates the case III with  = 10 when  = 10 −5 .Similarly, we illustrate the results of Theorem 6.1 using  = 0.75.The corresponding results are displayed in Figures5a-5c.The blue points in each picture in Figures4a-4c and 5a-5cshow the (sample) mean square approximation errors of N = 100 realizations of  , L().The red line in each figure shows the corresponding theoretical approximation upper bound.The numerical results show that the convergence rates of the mean  2 -error of   are consistent with the corresponding theoretical results in Theorem 6.1.