Homogenization of Fractional Kinetic Equations with Random Initial Data

We present the small-scale limits for the homogenization of a class of spatial-temporal random ﬁelds; the ﬁeld arises from the solution of a certain fractional kinetic equation and also from that of a related two-equation system, subject to given random initial data. The space-fractional derivative of the equation is characterized by the composition of the inverses of the Riesz po-tential and the Bessel potential. We discuss the small-scale (the micro) limits, opposite to the well-studied large-scale limits, of such spatial-temporal random ﬁeld. Our scaling schemes involve both the Riesz and the Bessel parameters, and also involve the rescaling in the initial data; our results are completely new-type scaling limits for such random ﬁelds.


Introduction
The purpose of this paper is to present new-type spatial-temporal scaling limits for the homogenization theory associated with a certain fractional kinetic equation and also a related system of reaction-diffusion type formed from it: where µ i > 0, 0 < β ≤ 1, 0 < α ≤ 2, 0 ≤ γ, and t > 0, x ∈ R n . The parameter β denotes the time-fractional index, and α, γ denote the space-fractional indices, for which we refer as the Riesz parameter and the Bessel parameter respectively (see [24, V.1 and V.3]).
When (β, α, γ) = (1, 2, 0), the system (1) is reduced to a classical reaction-diffusion system. The time-fractional index β < 1 means sub-diffusive (super-diffusive in case β > 1, which we do not study in this paper; see the remarks in Section 4). The spatial-fractional Riesz index α means the jumps of the evolution, and Bessel index γ means the tempering of large jumps; see the now-classic book of Stein [24, V.1 and V.3] for precise mathematical explanations. To our knowledge, fractional kinetic equations of Riesz-Bessel type appear firstly in Anh and Leonenko [1,2]; subsequent works in this direction by the authors and collaborators can be seen in [3,4,5,13,15,16] and the references therein. We mention that fractional differential equations involving more-than-one fractional parameters could be natural mathematical objects to describe long-range dependence and/or intermittency; one can find data exhibiting such characteristics in quite several fields including finance, hydrology, telecommunications, and turbulence.
In this paper, we aim to study the small-scale limit of spatial-temporal random field arising from the solution of the fractional Riesz-Bessel equation and its associated system (1) with µ 1 = µ 2 = µ > 0, subject to given the random initial data u 0 and v 0 , of which are independent and each one has a certain long-range dependence in its random structure. We mainly focus on the single equation case, and then the derive the system case via a certain decoupling method used in [17,18]. The opposite large-scale limits have been explored intensively in the above cited literatures; a notable point in such large-scale limits as in [1,2] is that the Bessel parameter does not play its role; while in our small-scale limits, both the Riesz and the Bessel parameters play their roles, and moreover we need the rescaling on the initial data. Such small-scale spatial-temporal scaling limit is a completely new result in the homogenization theory of random fields, to our best knowledge. Furthermore, our study may show how the theory of Riesz potentials and Bessel potentials, as in the Chapter V of Stein [24], may appear significantly in the homogenization of random fields.
The underlying idea in this paper is motivated by those works in [1,2,5,15,16] and the references therein. Namely, we use the spectral representations to describe the sample field arising from the initial data, and the relations between Hermite polynomials and homogeneous chaos associated with the initial data, to get representations for the limit field in terms of multiple Itô-Wiener integrals. From limit theorems point-of-view, our results, though it is the small-scale rather than the wellstudied large-scale as in those in the above citations, still could be regarded in the realm of noncentral limit theorems for convolution type integrals, in which the papers [25,11] are pioneering; see also the monograph of Major [22] and survey papers in the special volume edited by Doukhan, Oppenheim and Taqqu [9] The paper is organized as follows. In the preliminary Section 2, we list some result for the subsequent needs; we give the explicit solution of the system (1), which includes automatically the single equation solution, and our initial data are assumed to be the subordinated Gaussian fields. In the main Section 3, we present main results, the small-scale limit theorems of the solution field of the single equation and the related system, (1), in which we adopt usual time-derivative and the fractional spatial-derivative characterized by the Riesz and the Bessel parameters. The corresponding large-scale limits and time-fractional extensions are described in Section 4, together with two remarks for perspectives. The proofs of main results are give in the final Section 5.

Preliminaries
To begin with, we rewrite the system (1), with µ 1 = µ 2 = µ > 0, α > 0, γ ≥ 0 and β = 1, in the matrix form as follows: subject to some initial conditions The Green function G(t, x; α, γ) associated with the operator ∂ t + µ(I − ∆) γ 2 (−∆) α 2 , for which the spatial part is a hybrid of the Laplace and the Bessel operators, is defined via the spatial Fourier transform as follows; see, [24,Chapter 5] where < ·, · > denotes the inner product on R n . We remark that, by (4) it has In order to get an explicit representation for the solution of (2), we impose the following assumption on the matrix B. Condition DM. Suppose the matrix [b i j ] 1≤i, j≤2 is diagonalizable, i.e., the matrix B can be written as where P is a real-valued non-degenerate eigenvector matrix associated with the matrix B, and where d j is the eigenvalue associated with the eigenvector (p 1, j , p 2, j ) T (here and henceforth, T denotes the transpose). Without loss of generality, we suppose that d et(P) = 1. We always write Next, we consider the randomness of the initial data. Let (Ω, , ) to be an underlying probability space, such that all random element appeared in this paper are measurable with respect to it. In this paper, we assume that the initial data are a kind of isotropic stationary fields, namely the subordinated Gaussian fields; it is described as follows: Condition SGRID. We assume each component of the initial data (3) w 0 (x) := (u 0 (x), v 0 (x)) = (η 1 (x), η 2 (x)), x ∈ R n , is a random field of the following form in which ζ 1 (x) and ζ 2 (x) are independent, mean-square continuous, homogeneous and isotropic Gaussian random fields, each is of mean zero and of variance 1, and for each the spectral measure F j (dλ) has the (spectral) density f j (λ), λ ∈ R n , j ∈ {1, 2}, respectively. Each f j (λ) is decreasing for |λ| > λ 0 for some λ 0 > 0 and continuous for all λ = 0. Moreover, we assume that h j (·), j ∈ {1, 2}, are real non-random Borel functions satisfy Under the Condition SGRID, we have the spectral representation for the sample paths of ζ j (x), j ∈ {1, 2}, as below: where W j (·), j = 1, 2, are two orthogonal standard Gaussian noise measures on R n . Moreover, under (9), we can consider the following Hermite expansions of h j (u) in the Hilbert space L 2 (R, p(u)du) in which the Hermite coefficients The Hermite rank of the function h j (·) is defined by We also need the following two important properties ( see, for example, Major [22], Corollary 5.5 and p. 30) : and in the above, the integral representation (14) is the ρ-fold Itô-Wiener integral; the integration notation means that it excludes the diagonal hyperplanes z i = ∓z j , i, j = 1, ..., ρ, i = j.
We now can state the following explicit form of the spatial-temporal random fields of the system, subject to the Conditions DM and SGRID. It contains the single-equation solution consequently, which can be seen in [2, Section 2].
The above display of the solution can be obtained by firstly taking the spatial Fourier transform on both sides of (2). The Condition DM allows us to decouple the system; then we make the inverse Fourier transform to express the solution as the convolution of the Green function and the initial data. Then we make use of the Green function G(t, y; α, γ) defined as (4). Finally we plug in the initial data as described in the Condition SGRID.
We also impose the following assumption which is related to the long-range-dependence of the underlying Gaussian fields ζ j (x), j ∈ {1, 2}; we refer to [9] for the notion and the literatures of long-range dependence. In the following and henceforth, the notation f (·) ∼ g(·) means that the ratio f (·)/g(·) tends to 1, as the indicated variable "·" tends to ∞ or tends to 0, according to the context. Condition LD. The Gaussian random fields ζ j (x), j ∈ {1, 2}, in the Condition SGRID, have their covariance functions to be regular varying at infinity in the sense that: the covariance function R ζ j is given by where L : (0, ∞) → (0, ∞) is a slowly varying function at infinity and is bounded on each finite interval; recall that L is said to be slowly varying at infinity if lim y→∞ [L(c y)/L( y)] = 1 uniformly for any Here, the relation " ∼ " in (15) means that lim |x|→∞ |x| κ j R ζ j (x)/L(|x|) = 1; the similar notation will be used in this section and also in Section 5.
Under the Condition LD, by a Tauberian theorem (see, for example, the book of Leonenko [14, p. 66]), the spectral density functions of the random fields ζ j (x), j ∈ {1, 2}, are regular varying near the origin as follows: where K(n, κ j ) denotes the Tauberian constant. We note that, for each natural number ρ ≥ 2, the power of the covariance function (R ζ j (x)) ρ itself is still the covariance function of some random fields, for which there exists the corresponding spectral density function ( f j ) * ρ (λ), which is the ρ-th convolution of f j (λ): Note that L ρ (|x|) is still a slowly varying function for any ρ. Thus, when the ρ satisfies 0 < ρκ j < n, we can apply the Tauberian theorem again to get for j ∈ {1, 2}. While, if ρκ j > n then the covariance function (R ζ j (x)) ρ belongs to the class L 1 (R n ); thus the corresponding spectral density function is everywhere continuous and satisfies where The displays (17), (18) and (19) will be used in the proofs in Section 5.

Main results: the small-scale limits
In this section, we present the main results of this paper, which concerns with the small-scale (or say the micro) limits of the homogenization of the spatial-temporal random field described in Section 2. It features that both the Riesz parameter α and the Bessel parameter γ play their roles in the scaling procedure, and we also need to rescale the initial data.
Firstly, we state our result for the small-scale limit of a single fractional kinetic equation. The equation is then of the form: To our knowledge, the homogenization presented below is completely new spatial-temporal scalings.
In the below, the notation imposed on ζ wants to mean that the variable of ζ is under the indicated dilation factor. Theorem 1. Let s := s(t, x; s 0 (·)), t > 0, x ∈ R n , be the mean-square solution of (20), which satisfies the Conditions SGRID and LD, with κ ∈ (0, n m ), where m denotes the Hermite rank of the non-random function h(·) on R, which has the Hermite coefficients Section 2). Then, for any fixed parameter χ > 0, (1) the covariance function of the rescaled random field (2) When → 0, the rescaled random field s (t, x), t > 0, x ∈ R n , converges to the limiting spatialtemporal random field s m (t, x), t > 0, x ∈ R n , in the finite dimensional distribution sense, and s m (t, x) is represented by the multiple-Wiener integral where · · · denotes a m-fold Wiener integral with respect to the complex Gaussian white noise W (·) on R n .
Remark. The requirement for rescaling the initial data could be figured intuitively as that, the small-scale is meant to "freeze down" both the time t and the space x, to ensure the overall renormalization we have to "heat up" the initial data, so that it can be achieved a non-degenerate (though singular) limiting field. The term "small-scale limits", opposite to large-scale (other naming such as micro v.s. macro, local v.s. asymptotic), is adapted from the recent studies on the multi-scaling behavior of fractional non-Gaussian random fields, notably in [6,7]. It has also been used in the study of the fractal nature of the intersection local time measure [23].
Next, we state the small-scale limit for the system as follows. It shows that there is a nontrivial combination limit, namely the case (3) below, caused by the imposed relation on various parameters.
(1) If m 2 κ 2 > m 1 κ 1 , then the finite-dimensional distributions of the rescaled random field converge weakly, as → 0, to the finite-dimensional distributions of the random field with W 1 (·) is a complex Gaussian white noise on R n ( i.e., (23) with m, κ and W replaced by m 1 , κ 1 and W 1 , respectively).
(2) If m 1 κ 1 > m 2 κ 2 , then the finite-dimensional distributions of the rescaled random field converge weakly, as → 0, to the finite-dimensional distributions of the random field and W 2 (·) is a complex Gaussian white noise on R n ( i.e., (23) with m, κ and W replaced by m 2 , κ 2 and W 2 , respectively).
To understand the stochastic structure of the limiting fields, we state, for instance, the following covariance result of (Y * * * 1 (t, x) Y * * * 2 (t, x)).

Proposition 2.
For each fixed t > 0, the limiting vector field Y * * * 1 (t, x) Y * * * 2 (t, x) in the case (3) of Theorem 2 is spatial-homogeneous and its covariance matrix has the following spectral representation Remark. In view of the singularity of the (diagonal) spectral matrix near the origin, we may conclude that, for limiting vector field in the case (3), the long-range-dependence only exists within each individual component. In Proposition 3 below, we will find that there is also long-range-dependence between the different components of the large-scale limiting fields.

Large-scalings
In this subsection, we state the large-scale (or say the macro) limits of the system, in which only the Riesz parameter α plays its role in the scaling scheme. The result is compared to the single-equation case in And and Leonenko [2, Theorems 2.2 and 2.3]; it shows various limit fields may happen because of the different relations of the various parameters. The proof can be proceeded by the decoupling method.
(1) If m 2 κ 2 > m 1 κ 1 and d 1 > d 2 , then the finite-dimensional distributions of the rescaled random field converge weakly, as → 0, to the finite-dimensional distributions of the random field with W 1 (·) is a complex Gaussian white noise on R n .
(2) If m 1 κ 1 > m 2 κ 2 and d 1 > d 2 , then the finite-dimensional distributions of the rescaled random field converge weakly, as → 0, to the finite-dimensional distributions of the random field and W 2 (·) is a complex Gaussian white noise on R n .
The solution of (30) can be obtained via the fractional (both in time and in space) procedure (see, for example, [20,21]); under the Condition DM, we express the solution as the convolution of the fractional (both in time and in space) Green function and the initial data, as follows.
with the fractional Green function G β (t, x; d j ) is defined via the transformation where E β (·) is the Mittag-Leffler function defined by (see, for example, [2] or [8, For the properties of the Mittag-Leffler function, we refer the classic book by Erdélyi et.al. [10] (pp. 206-212, in particular p. 206 (7) and p. 210 (21)).
The following results are time-fractional versions of those in the above; however, the sub-diffusivity brings some new feature. For the large-scalings of homogenization of the system, we need to take an additional scaling on the matrix B in the system (30) in order to compromise the effect of this sub-diffusivity upon the interaction between u and v; while the sub-diffusivity has no influence on the small-scalings and thus it is the same as that in Section 3. To make the situation clear, in the following we denote the vector solution by w(t, x; w 0 (·), B). We only state some partial assertions in the below. Large-scalings: under m 1 κ 1 , m 2 κ 2 < min{2α, n} the finite-dimensional distributions of the rescaled random field t > 0, x ∈ R n , converge weakly, as → 0, to the finite-dimensional distributions of the random field , and for j ∈ {1, 2} T (1) (t, x; d j ) is expressed by the following multiple-Wiener integral: Small-scalings: under m 1 κ 1 , m 2 κ 2 < min{2(α + γ), n} the finite-dimensional distributions of the rescaled random field converge weakly, as → 0, to the finite-dimensional distributions of the random field

Two remarks for future study
1 Superdiffusive case. The time-fractional index β < 1 indicates the sub-diffusivity, and it changes to be the super-diffusivity if we consider β > 1 (see Section 1). In [18], the time-fractional reaction-wave type system with random initial data are studied, in which the first-order initial timederivatives of u and v play the crucial role. To consider spatial-temporal fractional kinetic systems which is super-diffusive in time and Riesz-Bessel in space will be a task of tremendous calculations. We also mention that, for the classical, i.e. non-fractional, heat-type system with random initial condition, the solution vector-field and the scaling limit are expressed in terms of heat kernels; this more explicit and simpler case is treated in [17].
2 Relativistic diffusion equation. That role of the Bessel parameter could be somewhat intriguing. It plays no role in the large-scaling [1,2], in which only the Riesz parameter appears; the Bessel parameter does appear in our small-scaling Theorem 1, and it plays jointly with the Riesz parameter. Perhaps, an appropriate viewpoint on this intriguing situation is to consider the "physically correct" Bessel operator, that is, to consider the prominent relativistic diffusion equation and its fractional version. This is addressed in [19].

Proofs of main results
We concentrate on the proofs of the main results in Section 3, that is, Theorem 1 and Theorem 2. The proofs of those propositions mentioned in Section 4 are skipped and left to the readers.
In the following proofs, ⇒ denotes the convergence of random variables (or random families) in distributional sense, and d = denotes the equality of random variables (or random families) in distributional sense. Moreover, we also denote f (t, x; ) g(t, x; ) if there exists a constant c := c(t, x) > 0 such that c g(t, x; ) < f (t, x; ) < c −1 g(t, x; ) when → 0.
The following lemma will play an important role in the proof of our main results; its proof is a easy consequence of the Cramer-Wold argument for two-dimensional random vectors (see, for example, [14, p. 6.]). U (t, x) and V (t, x) are independent random fields on R n × R + and Q (t) is a non-random 2 × 2 matrix. If there exist two random fields U 0 and V 0 such that U (t, x) ⇒ U 0 (t, x) and V (t, x) ⇒ V 0 (t, x), respectively, and Q (t) converges to Q(t) in the usual sense when → 0, then the finite dimensional distributions of X (t, x), t > 0, x ∈ R n , converge to the finite dimensional distributions of X :

Lemma 1. Let
Proof of Theorem 1 (1) Firstly, for simplification, we set N 1 ( ) = χ mκ L m ( −χ ). By the Hermite expansion, we can rewrite s( t, where the summation is in L 2 (Ω) sense. Hence, in accordion to the definition (21) about the random field s (t, x), it can be rewritten as with By the solution form given in Section 2, with d 1 = d 2 = 0, we have For the bracket above, by substituting t → t and λ → − 1 where we have used the self-similar property for Gaussian random measure on R n in the last equality. Therefore, by the orthogonal property for the Gaussian white noise, we can get where f * ρ (·) is defined in (17).
(i) For ρ ∈ N with mκ ≤ ρκ < n and any δ > 0, by (39) and (18), with and, by choosing δ small enough and (18), where the asymptotic equivalence is guaranteed by the uniform convergence theorem for the slowly varying function (see, for example, Leonenko [14,Section 1.4]). We remark that the f , defined as a spectral density function, is bounded outside of zero (while the singularity at zero is from the long-range-dependence assumption, as employed in subsection 3.2); therefore its ρ-th convolution remains to have, at most, the singularity only at zero. The conclusion of (i): Apart from the term Cov(I m (t, x)I m (t , x )), where we have used the fact that {l ∈ N|mκ ≤ łκ < n} is a finite set and on this set lim (ii) For ρ ∈ N with ρκ > n , by (39), (19) and Finally, from the expansion (36) for the random field s (t, x) and combining the observations (41) and (42)  Cov(s (t, x)s (t , x )) = (C m (h)) 2 K(n, mκ) R n e i<x−x ,τ> e −µ(t+t )|τ| α+γ |τ| n−mκ dτ.
(2) From the above discussion, we may apply Chebyshev's inequality to obtain that: Therefore, in view of Slutsky's argument (for example, see again [14, p. 6.]) , we suffice to focus our attention on the term I m (t, x). In the following we will prove that I m (t, x) converges in the distributional sense to s m (t, x), which is defined in (23), for each fixed (t, x) ∈ R + × R n . By the definition of N 1 ( ) and replacing the letter ρ by m in (38), we can rewrite (38) as follows with which, when → 0, satisfies lim →0 M (λ) (16) = (K(n, κ)) Now, applying the isometric property of the multiple Wiener integrals to the difference of (43) and (23), we have (44), the assumption f (λ) on Condition SGRID, and for mκ < n, by the Riesz composition formula (see, for example, [4, Appendix A]). Finally, the assertion (2) follows from Slutsky's and the Cramer-Wold arguments.

Proof of Theorem 2(1)
From the expression of the solution, we have where Q(t; d 1 , d 2 ) is the matrix defined by (7), and U(t, x) and V (t, x) are convolutions of the Green function with the respective initial data. By (45), We make use of the assumption that m 2 κ 2 > m 1 κ 1 .
Firstly, by Theorem 1(2), we have where X (1) Therefore, we may apply Lemma 1 to those U (t, x), V (t, x) and Q (t) on the above to obtain that

Proof of Theorem 2(2)
The proof is proceeded as the case (1), with the roles of m 1 , κ 1 , m 2 , κ 2 are interchanged, and thus we skip it.

Proof of Proposition 2
x) X (2) m (t, x) Because the representation for the limiting fields X (1) m (t, x) and X (2) m (t, x) is the same as the limiting field s m (t, x), defined in (23), we can apply the result (22) to get R n e i<x−x ,τ> e −µ(t+t )|τ| α+γ |τ| n−mκ dτ.