Fractional Riesz-Bessel equations extend classical kinetic models by using a composite fractional operator that combines the non-local behaviour of the fractional Laplacian with the regularising effect of the Bessel potential operator. Such equations are widely used in the study of transport in disordered media, multifractal dynamics in financial systems, and high-frequency turbulence, see [24, 37]. The work [19] initiated the development of this area, while [38] and [42] established detailed mathematical frameworks for fractional operators and integrals that are widely used today. Earlier publications considered these equations with random initial conditions given by Gaussian random fields [8, 9]. For such models, the asymptotic behaviour of the solution is determined not only by the differential operators but also by the spectral structure of the initial condition.

The first probabilistic approaches to the heat equation with random initial conditions were introduced by De Fériet and Rosenblatt [21, 44]. In recent years, there has been considerable interest in equations of this type with general random initial conditions. Several publications have demonstrated how the structure of these initial conditions influences the resulting solutions and their properties. The dependence structure of the initial random field is typically characterised by its spectral density: singularities at zero frequency correspond to classical long-range dependence, whereas singularities at nonzero frequencies are associated with cyclic long-memory behaviour [5]. Later work extended this analysis to broad classes of random environments and non-homogeneous settings [15, 47]. For classical diffusion equations, scaling limits of random solutions were derived in [2, 30], while related methods were applied to Burgers-type equations in [27, 29]. These studies established general principles for investigating asymptotics of rescaled random solution fields.

More recent research has focused on fractional operators and long-range dependence. Gay and Heyde [23] showed that fractional Laplace operators generate random fields with long-range dependence. Further studies of stochastic heat equations with fractional Laplace operators were carried out in [7, 13]. The results by Anh and Leonenko [8] proved that when the initial spectrum has a singularity at the origin, suitably normalised solutions converge to non-Gaussian limits. This established a direct link between fractional operators, singular spectral densities, and limit theorems. The analysis was further developed by renormalisation and homogenization methods for fractional-in-time/space diffusion equations with random input [9, 11]. These results belong to the non-central limit theory for strongly dependent random fields, see [22, 46]. Related asymptotic problems have been studied for Airy and Korteweg–de Vries equations [16, 25]. Stochastic models on spheres and hyperbolic diffusion equations were analysed in [14, 18]. Further results for diffusion in expanding space–time were recently obtained in [17].

Contrary to the classical increasing domain asymptotic settings, multiscaling refers to asymptotic zoom-in regimes in which the normalisation and the space–time rescaling are used to study local behaviour. For stochastic fractional diffusion equations, these settings lead to multiscaling limit theorems for their solutions, as well as the classification of the limit behaviour depending on the initial conditions [4]. The analysis of fractional Riesz–Bessel equations with cyclic long-memory initial conditions in [32] demonstrated that such rescaled solutions converge to Gaussian random fields. Related scaling and homogenization questions have also been studied for time-fractional relativistic diffusion equations and for coupled parabolic and reaction–diffusion–wave systems [33–35]. These works demonstrated how the limit behaviour depends on the fractional operator and the initial spectrum.

Kernel smoothing is often introduced to obtain improved behaviour and non-degenerate limits by controlling high-frequency contributions [4]. In numerous applications, this approach functions as a filter applied to realisations of random processes or fields [6, 12]. For example, the work of [50] considered kernel smoothing methods for stochastic partial differential equations, where solutions to the Fokker–Planck equation were approximated using adaptive Gaussian kernels. The approach was combined with stochastic filtering to estimate probability distributions from noisy observations. This established a direct link between kernel smoothing, types of considered equations and properties of the limits.

The existing literature (see [3–5] and the references therein) provides multiscaling theorems for several classes of equations with long-range and cyclically dependent initial conditions, including fractional and higher-order equations. However, the results in [3, 4] demonstrate that multiscaling limits do not exist in all cases and suggest using kernel smoothing to obtain non-degenerate asymptotic behaviour. To the best of our knowledge, no existing work has considered fractional Riesz–Bessel equations in conjunction with kernel smoothing and multiscale asymptotic analysis under general spectral settings. This study considers both classical and cyclic long-memory initial conditions, extending earlier results on the scaling limits of the corresponding solution fields. It also investigates the smoothness properties and dependence structures of these limit fields.

The paper is organised as follows. Section 2 provides the main definitions, notations, and some basic assumptions used throughout the paper. Section 3 derives limit theorems for filtered random fields with random initial conditions, for both classical and cyclic long-range dependence cases. It also presents numerical examples with realisations of the limit fields and their covariance structures. Hölder continuity and long-range dependence properties of the obtained multiscaling limit fields are studied in Section 4. Section 5 concludes the paper with a brief discussion and directions for future research.

All numerical computations and plotting in this paper were performed using the software R (version 4.6.0). The corresponding R code is freely available in the folder “Research materials” from the website https://sites.google.com/site/olenkoandriy/.

This section provides the main notations and background material required in the following sections.

Standard notations such as c ‾ for the complex conjugate of c, B ( R ) for the σ-field of Borel sets on R, Γ ( · ) for the Gamma function, and 1 A ( · ) for the indicator function of a set A are used throughout. The symbol C denotes constants, which exact values are not important for this exposition and may vary, even within the same proof.

For s ∈ R, the Sobolev space H s ( R ) consists of all functions f ∈ L 2 ( R ) such that ‖ f ‖ H s 2 : = ∫ R ( 1 + | λ | 2 ) s | f ˆ ( λ ) | 2 d λ < + ∞.

Let u ( t , x ) be a real-valued function of the two arguments x ∈ R d and t > 0. The time derivative of order β ∈ ( 0 , 1 ] is defined by ∂ β u ( t , x ) ∂ t β : = ∂ u ∂ t ( t , x ) , if β = 1 , D t β u ( t , x ) , if β ∈ ( 0 , 1 ) , where D t β u ( t , x ) : = 1 Γ ( 1 − β ) ∂ ∂ t ∫ 0 t ( t − τ ) − β u ( τ , x ) d τ − u ( 0 , x ) t β , t > 0 , is the fractional derivative in the Caputo–Djrbashian sense [42, (2.138)].

The publications [3, 8] and [9] studied fractional Riesz–Bessel equations (FRBE) (1) ∂ β u ( t , x ) ∂ t β = − μ ( I − Δ ) γ / 2 ( − Δ ) α / 2 u ( t , x ) , t > 0 , subject to a random initial condition (2) u ( 0 , x ) = ξ ( x ) , where α ≥ 0, γ > 0, and μ > 0, while ξ ( x ), x ∈ R d , is a zero-mean Gaussian random field. For a function φ in the Schwartz space S ( R d ), with Fourier transform φ ˆ ( λ ) : = ∫ R d e − i ⟨ x , λ ⟩ φ ( x ) d x, the fractional Laplacian ( − Δ ) α / 2 is defined by ( − Δ ) α / 2 φ ˆ ( λ ) = | λ | α φ ˆ ( λ ), and the Bessel fractional differential operator ( I − Δ ) γ / 2 is defined by ( I − Δ ) γ / 2 φ ˆ ( λ ) = ( 1 + | λ | 2 ) γ / 2 φ ˆ ( λ ). Thus, ( − Δ ) α / 2 is the inverse of the Riesz potential operator ( − Δ ) − α / 2 , whose Fourier multiplier is | λ | − α , while ( I − Δ ) γ / 2 is the inverse of the Bessel potential operator ( I − Δ ) − γ / 2 , whose Fourier multiplier is ( 1 + | λ | 2 ) − γ / 2 . When α = 0, the operator ( − Δ ) α / 2 is understood as the identity operator. See, for example, [8, 9].

The Cauchy problem (1) with the non-random initial condition given by the Dirac delta function, u ( 0 , x ) = δ ( x ), has a unique solution determined by G ˆ ( t , λ ), t > 0, λ ∈ R d , the Fourier transform with respect to the second argument of the fundamental solution (i.e., the Green function) of the non-random Cauchy problem.

The solution to the initial value problem (1) and (2) in convolution form is defined as (3) u ( t , x ) = ∫ R d G ( t , x − y ) u ( 0 , y ) , d y , where the Green function G ( t , x ) is specified by its Fourier transform.

In the literature this solution is called a Green-type solution or a mean-square solution, as it can be interpreted in the mean-square sense. For more details and justifications of this approach, consult [10, 12].

For simplicity, this paper focuses on the one-dimensional case of d = 1. To make the notations consistent with multidimensional cases considered in the mentioned publications, the operator Δ = ∂ 2 / ∂ x 2 denotes the one-dimensional Laplacian.

Let ξ ( x ), x ∈ R, be a measurable stationary real-valued zero-mean Gaussian random process defined on a probability space ( Ω , F , P ). Let Z ( · ) and W ( · ) denote respectively, an orthogonal random measure of ξ ( · ) and a Gaussian white-noise random measure on R, with E | Z ( d λ ) | 2 = F ( d λ ), where F ( · ) is the spectral measure of ξ ( · ). In this paper, the random measures Z ( · ) and W ( · ) are assumed to be symmetric in the sense of complex conjugate symmetry, that is, Z ( A ) = Z ( − A ) ‾ and W ( A ) = W ( − A ) ‾ for all A ∈ B ( R ), which ensures that the process ξ ( x ) is real-valued.

The covariance function of ξ ( x ) can be written as (4) r ( x ) : = Cov ( ξ ( 0 ) , ξ ( x ) ) = ∫ R e i x λ F ( d λ ) , where F ( · ) is the spectral measure.

If the spectral measure is absolutely continuous, it can be represented as F ( Δ ) = ∫ Δ f ξ ( λ ) d λ , Δ ∈ B ( R ) , where the function f ξ ( λ ) , λ ∈ R is called the spectral density function of the stationary process ξ ( · ). Since the stationary process ξ ( · ) is real-valued, its spectral density is necessarily even.

Then, the following spectral representation of the random process ξ ( x ) holds true ξ ( x ) = ∫ R e i λ x Z ( d λ ) = ∫ R e i λ x f ξ ( λ ) W ( d λ ) .

The Fourier transform of the Green function G ( t , x ) is equal to, see [10, 12], G ˆ ( t , λ ) = E β − μ t β | λ | α 1 + | λ | 2 γ / 2 , Here E β ( · ) denotes the Mittag–Leffler function, defined by E β ( s ) : = ∑ k = 0 ∞ s k Γ ( 1 + β k ) , s ∈ R , 0 < β ≤ 1 . For β = 1, E 1 ( s ) = e s . More details on the Mittag–Leffler function and its properties can be found in [36].

For 0 < β < 1 and negative values of its argument, the Mittag–Leffler function satisfies the following two-sided estimate (see [45, Theorem 4, inequality (6.8)]): (5) 1 1 + Γ ( 1 − β ) s ≤ E β ( − s ) ≤ 1 1 + s Γ ( 1 + β ) , s ≥ 0 .

The solution (3) of the initial value problem (1) and (2) admits the following stochastic integral representation, see [12], (6) u ( t , x ) = ∫ R e i x λ E β ( − μ | λ | α ( 1 + | λ | 2 ) γ / 2 t β ) Z ( d λ ) .

The covariance function of the solution field u ( t , x ) is Cov ( u ( t , x ) , u ( t ′ , x ′ ) ) = ∫ R e i ( x − x ′ ) λ E β ( − μ | λ | α ( 1 + | λ | 2 ) γ / 2 t β ) × E β ( − μ | λ | α ( 1 + | λ | 2 ) γ / 2 ( t ′ ) β ) F ( d λ ) .

To further specify the class of Gaussian random processes ξ ( x ), x ∈ R, used as the initial condition (2), the following assumption is used, see [5] for more details.

The frequency parameters w j determine the locations of the singularities of the spectral density. In particular, for each j = 1 , … , n with A j > 0, the corresponding singularities are located at λ = ± w j , while, when A 0 > 0, w 0 = 0 corresponds to a singularity at the origin.

Let us define the next constants C ( κ ¯ , w ¯ , A ¯ ) : = ∑ j = 1 n c 1 , j ( κ j ) A j K κ j − 1 2 ( | w j | ) | w j | κ j − 1 2 , C ( κ 0 , A 0 ) : = A 0 2 Γ κ 0 cos ( κ 0 π / 2 ) , where c 1 , j ( κ j ) : = 2 1 { 0 } ( j ) 2 1 − κ j 2 / π Γ κ j / 2 , j = 0 , 1 , … , n, κ ¯ = ( κ 1 , … , κ n ), w ¯ = ( w 1 , … , w n ), and A ¯ = ( A 1 , … , A n ). The symbol K ν ( · ) denotes the modified Bessel function of the second kind, which can be defined, see [49, p. 78], by K ν ( z ) : = 1 2 ∫ 0 + ∞ s ν − 1 exp − 1 2 s + 1 s z d s , z > 0 , ν ∈ R .

The covariance function in (7) is non-integrable and exhibits an oscillating behaviour, which corresponds to the cyclic long-range dependence scenario. It follows from (7) that the corresponding spectral density has the representation (8) f ξ ( λ ) : = ∑ j = 1 n c 1 , j ( κ j ) 2 A j ( K κ j − 1 2 | λ + w j | | λ + w j | κ j − 1 2 + K κ j − 1 2 | λ − w j | | λ − w j | κ j − 1 2 ) + c 1 , 0 ( κ 0 ) 2 A 0 K κ 0 − 1 2 ( | λ | ) | λ | κ 0 − 1 2 ,

The publications [4, 5] showed that, for the Cauchy problem with random initial conditions, there are cases where multiscaling limits do not exist. They also proposed kernel smoothing as a way to obtain non-degenerate limits in the case of higher-order heat equations. The following results apply this approach to the fractional Riesz–Bessel equations studied in [3]. For the same class of random initial conditions, the multiscaling limit results in [3] were obtained under the restriction α > 1 / 2 when A 0 = 0 and under the restriction α > κ 0 / 2 when A 0 > 0. The filtering approach developed in the present paper provides additional integrability through the kernel multiplier and allows these restrictions on α to be removed.

Consider a spatial kernel g ( x 1 , x ), x 1 , x ∈ R, and define the centered spatially averaged random field (9) U ε g ( t , x ) : = ε ρ 1 ∫ R g ( x 1 , x ) u t ε ρ 2 , x 1 ε ρ 3 d x 1 , where ρ i ∈ R, i = 1 , 2 , 3, are some scaling parameters. As the initial condition is given by a zero-mean Gaussian random field, the fields u ( t , x ) and U ε g ( t , x ) are also zero-mean.

We will use real-valued spatial kernel functions g ( x 1 , x ), x 1 , x ∈ R, such that for all x it holds g ( · , x ) ∈ L 1 ( R ) ∩ L 2 ( R ). Their Fourier transforms with respect to the first argument are given by g ˆ ( λ , x ) : = ∫ R e − i λ x 1 g ( x 1 , x ) d x 1 . To obtain non-degenerate limits, we assume that g ˆ ( · , x ) ≠ 0 on a set of positive Lebesgue measure. Since the Fourier transform of these kernel functions is taken with respect to the first variable, we will use the notation ‖ g ( z , x ) ‖ H z s for the corresponding Sobolev norm in the z-variable. Theorem 1.

Consider the random field u ( t , x ), t > 0, x ∈ R, defined by (1) and the random initial condition given by (2), satisfying Assumption 1 with A 0 = 0.

For ρ 1 = − ρ 3 / 2, ρ 2 = α ρ 3 / β, and ρ 3 > 0, when ε → 0, the finite-dimensional distributions of U ε g ( t , x ) converge weakly to those of a zero-mean Gaussian random field  U 0 g ( t , x ) given by (10) U 0 g ( t , x ) : = C ( κ ¯ , w ¯ , A ¯ ) ∫ R g ˆ ( λ , x ) E β − μ t β | λ | α W ( d λ ) .

The limit random field U 0 g ( t , x ) has the covariance function (11) Cov ( U 0 g ( t , x ) , U 0 g ( t ′ , x ′ ) ) = C ( κ ¯ , w ¯ , A ¯ ) ∫ R g ˆ ( λ , x ) g ˆ ( λ , x ′ ) ‾ × E β ( − μ t β | λ | α ) E β ( − μ ( t ′ ) β | λ | α ) d λ .