Fractional Diffusion to a Cantor Set in 2D

: A random walk on a two dimensional square in R 2 space with a hidden absorbing fractal set F µ is considered. This search-like problem is treated in the framework of a diffusion–reaction equation, when an absorbing term is included inside a Fokker–Planck equation as a reaction term. This macroscopic approach for the 2D transport in the R 2 space corresponds to the comb geometry, when the random walk consists of 1D movements in the x and y directions, respectively, as a direct-Cartesian product of the 1D movements. The main value in task is the ﬁrst arrival time distribution (FATD) to sink points of the fractal set, where travelling particles are absorbed. Analytical expression for the FATD is obtained in the subdiffusive regime for both the fractal set of sinks and for a single sink.


Introduction
In this paper, we study a 2D continuous subdiffusive transport in the presence of an absorbing fractal set. This macroscopic consideration of the problem is treated in the framework of a diffusion-reaction equation, when an absorbing term is included inside a Fokker-Planck equation as a reaction term, and this approach is analogous to 1D search problems [1][2][3][4]. In the 2D case, this macroscopic continuous search task is less studied, although a theory of intermittent search tasks in 2D is well developed for both lattices [5][6][7] and continuous spaces [8,9].
The 2D movements in the xy-space consist of a direct product of 1D movements in the x and y directions separately. That is, the movements are topologically constrained to a comb geometry, while, at the same time, any point on R 2 is accessible due to this "search" strategy. The latter is a 2D random walk in the presence of absorbing sinks that eventually is reduced to the non-Markovian transport in a comb geometry (see Figure 1). A comb model is a simple caricature of various types of natural branched structures, where random walks on comb structures provide a geometrical explanation of anomalous diffusion. An extended review of a variety of realizations of comb models can be found in [10]. As a specific motivation for the present study, this model can be applied in tubular search processes like neurons [5,[11][12][13][14][15][16][17]. In particular, an intermittent search process on a dendrite tree has been considered [11], while a relation between the comb model and a spine dendrite structure has been established as well [12][13][14][15] with further explorations of anomalous diffusion inside spine dendrites [16][17][18]. Another important application relates to fabrication of fractal absorbers in metamaterials applied for antenna technology (see, e.g., [19] and references therein). A random walk on a comb in the presence of a fractal sink structure. This fractal structure, where the sink points are placed at y = Y, is shown in zoom. The comb consists of the backbone-x axis and the y direction corresponds to continuously distributed fingers along the backbone. This structure results from the phenomenological consideration of 2D diffusion with diffusion coefficients D xx = D x δ(y) and D yy = D y , where D x , D y are constant values [20]. The main question in task is the first arrival time distribution (FATD) for a fractal structure of sinks F µ , where traveling particles are absorbed. A general analytical expression has been obtained, and we show that, in the subdiffusive regime, the FATD has an universal time behavior ∼t −3/4 , when the power law index is independent of the fractal dimension µ in the large time asymptotic. The suggested analytical method makes it possible to study the FATD for anomalous diffusion in the presence of a single sink, as well. It is well known that the 2D diffusion in R 2 is recurrent [21]. Therefore, an analytical expression for the FATD for a single sink, say at a point (x, y) = (X, Y), is another important question in task. However, a straightforward solution of the macroscopic equation with µ = 0 cannot be obtained, in contrast to the case with µ = 0, which can be solved analytically. Therefore, first, we use the analytical solution for the FATD for a fractal set of sinks with µ = 0. Then, taking into account that the fractal set in continuous space is not a frozen-fixed structure, the limit µ → 0 is performed, and the FATD for a single sink is obtained in the form of the Fox H-functions.

Preliminaries: Subdiffusion and Boundary Conditions on a Comb
Important parts of the 2D continuous movements are boundaries and boundary conditions [8]. Therefore, in this section, we briefly outline some properties of the comb model, related to the finite boundary conditions. In Appendix B, we present more details of the fractional comb motion, when the boundaries are at infinity. The necessary elements of fractional calculus for this issue are presented in Appendix A. A macroscopic consideration of diffusion on comb structures is based on a Fokker-Planck equation [20]. According to the phenomenological approach of Arkhincheev and Baskin [20], the comb geometry implies that the displacement in the x-direction is possible only along the structure axis, i.e., along the x-axis at y = 0. Diffusion in the x-direction is highly inhomogeneous, and the diffusion coefficient is given by means of the Dirac delta function, D xx = D x δ(y), while the diffusion coefficient in the y-direction, the side-branch direction, also called teeth or fingers, is a constant, namely D yy = D y . This situation also implies that fingers are continuously and homogeneously distributed along the backbone. Then, the diffusion equation on the comb structure reads [20] ∂ ∂t P(x, y, t) = D x δ(y) ∂ 2 ∂x 2 P(x, y, t) + D y where P(x, y, t) is the probability density function (PDF) of finding a diffusing particle at time t at the position with coordinates (x, y) in the 2D comb space. The initial condition is P(x, y, t = 0) = δ(x)δ(y) and natural boundary conditions are chosen for the backbone at infinity, i.e., P(x, y, t) = ∂ x P(x, y, t) = 0 for x = ±∞, and reflecting boundary conditions are at y = ±h, namely ∂ y P(x, y, t)| y=±h = 0, for fingers.
Following standard procedures of the Fourier transformation with respect to the space variables and the Laplace transformation with respect to the time, LF x P(x, y, t)] =P(k, y, s), particular solutions of the comb equation can be obtained. In the Fourier-Laplace space, Equation (1) reads where s is the Laplace variable and k is the Fourier variable. We look for solutions in the form [22] P(k, y, s) = g 1 (k, s)g 2 (y, where g 1 (k, s) =P(k, y = 0, s) describes diffusion along the backbone in Laplace space. This also supposes that g 2 (y = 0, s) = 1 is the second boundary condition for fingers. In sequel we consider subdiffusion along the fingers, which takes place at the condition h s/D y 1 [22]. Substituting solution (3) into Equation (2) and integrating over y, we obtain Analytical solutions of the distribution functions can be represented in the form of the Mittag-Leffler and the Fox H functions, as described in Appendix B. It turns out that the transport along the backbone is subdiffusive until time t < T = h 2 /D y , when the boundaries are not reached by particles during normal diffusion, which occurs along the side-branched fingers.
For sequel analysis, it is convenient to work with dimensionless space-time variables and parameters. It can be done by a scaling procedure as follows: t · where the diffusion coefficients in both x and y directions are equaled to one. The dimensionless subdiffusive time scale is T = h 2 D 2 y /D 2 x .

First Arrival Time Distribution for Subdiffusion
As admitted above, the random walk in the presence of sinks is analogous to a search problem. An important characteristic of the search process is a first arrival time distribution (FATD) ℘ fa (t), which is an important characteristics of first-passage processes [21,23]. In a random search equation, it plays a role of a strength of a δ-sink term. For example, a Brownian random search process is governed by the diffusion equation with the δ-sink term in the form of the Dirac delta function with an unknown strength ℘ fa (t) at point (x, y) = (X, Y). Then, amending dimensionless Equation (5) with this δ-sink term of strength ℘ fa (t), we obtain the diffusion-reaction equation for the non-normalized PDF f (x, y, t), The initial condition is P(x, y, t = 0) = δ(x − x 0 )δ(y − y 0 ). Considering subdiffusion along the backbone, natural (zero) boundary conditions are chosen at infinity for both the backbone and fingers that ensures a subdiffusive transport along the backbone.
Integration over the R 2 space yields that the FATD relates to the survival probability S(t) as follows The random walk is finished when a particle arrives at the point (X, Y), where it is absorbed by the sink term, and the PDF vanishes, f (X, Y, t) = 0. Then, the condition f (X, Y, t) = 0 leads to the equation for the FATD ℘ fa (t). This condition is also valid in the Laplace spacef (X, Y, s) = 0.

Fractal Set of Sinks
As admitted above, our aim is to to find the time distribution of first arrivals at the sink term at point (X, Y). However, as shown in the sequel, a single point sink cannot be found in the 2D search, in contrast to the 1D search (when Y = 0). Therefore, to overcome this deficiency of the macroscopic study of the FATD, we introduce a fractal set of sinks as follows.
Let in the vicinity of an unknown (hidden) point (X, Y) there be a fractal set of sink points. For example, it is a random Cantor set with the fractal dimension µ The set is embedded in an -segment/interval R 1 ( ), that is F µ ⊆ R 1 ( ), where an arriving particle is absorbed by sinks.
Therefore, a particle with the initial point (x 0 , y 0 ) is absorbed arriving at any point where for simplicity of the notation ∑ j ≡ ∑ X j ∈F µ is used for sum over F µ (x) . As admitted above and is understood in the sequel analysis, the FATD ℘ fa (t) can be found only for a set of x points with nonzero (at least) fractal measure (Hausdorff, or box-counting) for the present geometry. Note that, otherwise, the equation in task becomes singular due to the last term (without the sum), and, as a result of this deficiency, the FATD cannot be found. Performing the Fourier transforms with respect to x and y, we obtain This yieldsf To obtain this equation in a closed form, we expressf (k x , y = 0, s) in terms of the known components off (k x , k y , s). To this end, we perform the Fourier inversion of Equation (11) with respect to y. This procedure yieldŝ and for y = 0 we obtain the relation Substituting Equation (13) into Equation (12) and performing the Fourier inversion with respect to k x , we obtain The obtained result readŝ Then, arriving at the point (x, y) = (X j , Y), a particle is absorbed andf (X j , Y, s) = 0 in Equation (15). First, however, the summation over the fractal set F µ is performed. One should recognize that fractal sets (like a Cantor set) are uncountable. Therefore, the sum ∑ j . . . in Equation (15) is purely formal and its mathematical realization corresponds to integration to fractal measure M µ ∼ l µ such that ∑ l i ∈F µ δ(l − l i ) = 1 Γ(µ) l µ−1 is the fractal density [24][25][26], and dM µ = 1 Γ(µ) l µ−1 dl. Here, Γ(µ) = (µ − 1)! is the gamma function. Therefore, summation over all points of the fractal set corresponds to integration with fractal measure, ∑ X i ∈F µ g(X i ) = 0 g(l)dM µ (l). Accounting for all i and x = X j that X i = X + l i and X j = X + l, where l, l i ∈ (0, ), we obtain Here, γ(µ, z) = Γ(µ) − Γ(µ, z) is an incomplete gamma function. For z 1, it is estimated by Γ(µ)z µ [27].
Taking into account Equations (16a) and (16b) and the sink condition f (X, Y, s) = 0, one obtains from Equation (15) In this limit, the second term in the denominator can be neglected, while the first one is of the order of ∼ √ s(|Y| + |y 0 |). In the nominator, only the first, largest exponential is taken into account. The Laplace inversion of Equation (18) yields the FATD for the fractal set F µ in the large time asymptotics. To that end, we present the exponential in the form of the Fox H function [28,29] which eventually yields . (20) Equation (20) is the main result, which determines the FATD of anomalously diffusive particles to the fractal set of sinks. Since asymptotic expansion of the Fox H function at zero is O(1), the time asymptotic behavior of the FATD is ℘ fa (t) t − 3 4 . The factor 1−µ Γ(µ) is a part of a stationary probability, and it yields the correct unite limit for µ = 1. The limit of a single sink point for µ → 0 does not exists in solution (20). However, 2D diffusion is recurrent, therefore a single sink can be found as well. This means that the limits µ → 0 and s → 0 do not commute. Then, to obtain the FATD for a single sink, the limit µ → 0 should be taken first in Equation (17).

Single Sink
Now, the limit is µ → 0, which means that the only one sink point is considered first, and then the limit of the large times s → 0 are considered in Equation (17). Eventually, for µ = 0, the absorbing point (X j , Y) = (X, Y) yields the Laplace image of the FATD,℘ fa (s) as followŝ 2|X−x 0 |s where ∆ y = |y 0 | − |Y| > 0. We believe that this restriction results from a correlation between the initial condition and the sink term. This restriction is discussed in Section 3. Performing the inverse Laplace transform of Equation (21), we present the result in the convolution form as follows That is, the FATD consists of two physically justified parts. The first one G y (t) is the FATD for the Brownian search along the fingers, when the worker passes the distance |y 0 | from the initial point at y 0 to the backbone at y = 0, while the second Brownian distance |Y| is a random walk from the backbone to the cantor set F ν (y). As expected, this Brownian part of the FATD corresponds to the well known Lévy-Smirnov distribution [30,31], The second part of the FATD, G x (t), corresponds to subdiffusion along the backbone [4]. Therefore, to estimate the FATD by the Laplace inversion, we follow the same procedure of Equations (19) and (20) that yields

Discussion
In the paper, we are concerned with a random walk on a 2D square S 2 (h × h) ∈ R 2 with a hidden absorbing fractal set F µ ⊆ R 1 ( ) ⊂ R 1 x and reflecting boundaries. This search-like problem is treated in the framework of a diffusion-reaction equation with the absorption term, which is included inside the Fokker-Planck equation as the reaction term. This macroscopic approach to the 2D movements in the R 2 space is considered in the comb geometry, when the random walk consists of the 1D movements in the x and y directions, respectively, as a direct-Cartesian-product of the 1D movements, where diffusion in fingers affects strongly the backbone transport. The main value in task is the first arrival time distribution (FATD) ℘ fa (t) of particles traveling to sink points, where they are absorbed. The imposed finite boundary conditions lead to the two time scales. The first, initial time scale, T with t < T = D 2 y h 2 /D 2 x , corresponds to subdiffusive regime, when the boundary conditions can be considered at infinity. In this case, the main result for the FATD for the fractal set of sinks in Equation (20) is not restricted, since t < T = ∞. The FATD for a single sink (22) is restricted according to the condition ∆ y = |y 0 | − |Y| > 0 and therefore by t < T = T(h). We however consider h as large as necessary to have a good statistics for the FATD. Moreover, taking for the initial condition y 0 = h, the sink point can always be found, since the condition for the inverse Laplace transform is always fulfilled due to the condition |y 0 | − |Y| = h − |Y| > 0. The restriction ∆ y = |y 0 | − |Y| > 0 results from subdiffuion. The latter results from the condition h → ∞, which leads for the first time scale t < T → ∞, and the time restriction for the FATD is lifted. Correspondingly, the sink point is always reachable, and the escape probability (the probability to be absorbed) is E = lim T→∞℘fa s = 1 T = lim T→∞ T 0 ℘ fa (t)dt = 1. As conclusive remarks, we also come to grips with some details of the escape probability E for the second time scale, t > T. It can be estimated from the reflecting boundary conditions of Equation (1), when both backbone and fingers are taken into account. Then, the Laplace image obtained from Equation (3) readŝ where the main contribution is due to the term with k = 0 [32]. The only outgoing flux is due to the sink. According to the comb geometry, the current has the y component j y (X, Y, t) = −∂ y P(x = X, y, t)| y=Y , only. Therefore, the hitting probability [21] and the FATD are determined by the current at (x, y) = (X, Y) and these values are The question on the FATD in the h × h square is open for the continuous macroscopic consideration. Nonetheless, it also follows from Equations (25), (26a) and (26b) that subdiffusion in the bounded backbone and opened fingers [16] is an efficient way to arrive at the sink. In this case when h → ∞, the solutions (25) and (26) yield (not tempered) the Lévy-Smirnov distribution (23) for the first passage probability and the escape probability is E =℘ fa (s = 0) = 1. Since the current is a continuous value it can be associated with the FATD ℘ fa (t) as well.
Author Contributions: A.I. and T.S. contributed to each part of this work equally. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.
where α > 0, x > a and Γ(z) is the Gamma function. Fractional derivation was developed as a generalization of integer order derivatives and is defined as the inverse operation to the fractional integral. Therefore, the fractional derivative is defined as the inverse operator to I α a x , namely D α a x f (x) = I −α a x f (x) and I α a x = D −α a x . Its explicit form is This integral diverges for arbitrary α > 0, and a regularization procedure is introduced with two alternative definitions of D α a x . For an integer n defined as n − 1 < α < n, we obtain the Riemann-Liouville fractional derivative, There is no constraint on the lower limit a. For example, when a = 0, one has This fractional derivation with the fixed low limit is also called the left fractional derivative. However, one can introduce the right fractional derivative, where the upper limit a is fixed and a > x. For example, the right fractional integral is Another important property is D α I β = I β−α , where subscripts are omitted for brevity. Note that the inverse combination is not valid. In general, I β D α = I β−α , since it depends on the lower limits of the integrations [33]. We also use here a convolution rule for the Laplace transform for 0 < α < 1 Note that for arbitrary α > 1 the treatment of the Caputo fractional derivative by the Laplace transformation is more convenient than that of the Riemann-Liouville one.
Note that the solutions considered here can be obtained in the form of the Mittag-Leffler function by the Laplace inversion [28,33,35] where C is a suitable contour of integration, starting and finishing at −∞ and encompassing a circle |s| ≤ |z| 1 ν in the positive direction, and ν, β > 0.

Appendix B. Fractional Fokker-Planck Equation
Integration with respect to y leads to a equation for the marginal distribution P 1 (x, t) = ∞ −∞ P(x, y, t)dy: whereP(x, y = 0, s) = f (x, s). To obtain this equation in closed form, we return to Equation (3) in (x, y)-space,P (x, y, s) =P(x, y = 0, s) exp − s/D y |y| , and integrate it with respect to y. This yields the relationP(x, y = 0, s) = s/4D yP1 (x, s). Using this relation and Laplace transforming Equation (A8), we find Performing the Laplace inversion, we obtain the integro-differential equation (A11) The integro-differential operator D C α 0 t ≡ D C α t is the so-called Caputo fractional derivative, while Equation (A11) is called the fractional Fokker-Planck equation (FFPE). Instead of the Caputo derivative, it is possible to employ the Riemann-Liouville fractional derivative, which leads to a different form of the FFPE [36], where α = 1/2. For arbitrary 0 < α < 1, this equation is a general form of the FFPE with the solution (A24) 2) .
In that case, the diffusion coefficient must be generalized as well, D x /(2 D y ) = D 1 2 → D α .
Using property of Equation (A22), we obtain