Entropy Production Rates of the Multi-Dimensional Fractional Diffusion Processes

Our starting point is the n-dimensional time-space-fractional partial differential equation (PDE) with the Caputo time-fractional derivative of order β,0<β<2 and the fractional spatial derivative (fractional Laplacian) of order α,0<α≤2. For this equation, we first derive some integral representations of the fundamental solution and then discuss its important properties including scaling invariants and non-negativity. The time-space-fractional PDE governs a fractional diffusion process if and only if its fundamental solution is non-negative and can be interpreted as a spatial probability density function evolving in time. These conditions are satisfied for an arbitrary dimension n∈N if 0<β≤1,0<α≤2 and additionally for 1<β≤α≤2 in the one-dimensional case. In all these cases, we derive the explicit formulas for the Shannon entropy and for the entropy production rate of a fractional diffusion process governed by the corresponding time-space-fractional PDE. The entropy production rate depends on the orders β and α of the time and spatial derivatives and on the space dimension n and is given by the expression βnαt, t being the time variable. Even if it is an increasing function in β, one cannot speak about any entropy production paradoxes related to these processes (as stated in some publications) because the time-space-fractional PDE governs a fractional diffusion process in all dimensions only under the condition 0<β≤1, i.e., only the slow and the conventional diffusion can be described by this equation.


Introduction
One of the most prominent and broadly-recognized applications of Fractional Calculus (FC) is for description of the anomalous transport processes [1][2][3][4]. The basic FC model for the slow anomalous diffusion is the time-fractional partial differential equation (PDE) that interpolates between the time-independent Poisson equation and the diffusion equation. In Reference [5], the fundamental solution to the time-fractional diffusion equation with the Caputo fractional derivative of order β ∈ (0, 1) and the spatial Laplace operator was shown to be non-negative and normalized. Thus it can be interpreted as a spatial probability density function evolving in time that provides a strong justification for the time-fractional diffusion equation with the time derivative of order β ∈ (0, 1) to act as a model for a diffusion process. This process is anomalous and slow because the mean squared displacement of the diffusing particle behaves as c β t β , β ∈ (0, 1) in contrast to the linear dependence of the mean squared displacement on time in the case of the conventional diffusion.
In analogy to the case of the slow anomalous diffusion, some FC models for the supper-diffusion (fast diffusion) processes were introduced in form of the time-, space, or time-space-fractional PDEs. In particular, the time-fractional diffusion-wave PDE with the Caputo fractional derivative of order β ∈ (1, 2] and the spatial Laplace operator, which interpolates between the diffusion equation and the wave equation, was often discussed in the literature in connection with the supper-diffusion. However, it turned out that the fundamental solution to the n-dimensional time-space-fractional PDE with the Caputo time-fractional derivative of the order β, 0 < β < 2 and the fractional spatial derivative (fractional Laplacian) of the order α, 0 < α ≤ 2 fails to be non-negative for n ≥ 2 for all 1 < β < 2, 0 < α ≤ 2 [6]. This means that for n ≥ 2 this equation with the time-fractional derivative of order β greater than one cannot serve as a model for any diffusion process but rather as a model for the damped waves propagation [7,8].
However, the situation in the one-dimensional case (n = 1) is very different from the one in the multi-dimensional case (n ≥ 2). As shown in Reference [9], the fundamental solution to the one-dimensional time-space-fractional PDE with the time-fractional Caputo derivative of order β ∈ (0, 2] and the fractional spatial Riesz-Feller derivative of order α ∈ (0, 2] and skewness θ (|θ| ≤ min {α, 2 − α}) is non-negative both for 0 < β ≤ 1, 0 < α ≤ 2 and for 1 < β ≤ α ≤ 2, i.e., this equation can describe some diffusion processes also for β greater than one. The entropy and the entropy production rates of the processes governed by the one-dimensional time-fractional diffusion equation with the Caputo time-fractional derivative of order β ∈ [1, 2] and the second spatial derivative (Laplace operator in the one-dimensional case) were discussed in [10,11]. According to References [10,11], the entropy production rate of the time-fractional diffusion equation increases with increasing of β from one (diffusion) to two (wave propagation) that results in the so called entropy production paradox. Many efforts were put into "resolving" of this paradox (see e.g., [12] and references therein). However, as mentioned above, this paradox does not appear in the dimensions two and three, i.e., in the cases that are important for applications. The paradox in the one-dimensional case is rather a mathematical caprice than a physical phenomenon and thus it has only very restricted-if any-relevance for applications.
In the literature, the entropy and the entropy production rates of the processes governed by some other particular cases of the time-space-fractional PDE with the Caputo time-fractional derivative of order β, 0 < β < 2 and the fractional spatial derivative (fractional Laplacian) of order α, 0 < α ≤ 2 have been considered. In References [13,14], the case of the one-dimensional space-fractional PDE with the first-order time derivative and fractional spatial derivative of order α, α ∈ (0, 2] was analyzed in detail. In Reference [15], a closed form formula for the Shannon entropy of the fundamental solution to the one-dimensional neutral-fractional PDE with the fractional derivatives of the same order α, 1 ≤ α ≤ 2 both in space and in time, was derived. The n-dimensional case of this equation was analyzed in Reference [7]. In Reference [16], it was shown that the entropy production rate of the fundamental solution to the one-dimensional α-fractional diffusion equation is exactly the same as in the case of the conventional diffusion equation. The α-fractional diffusion equation is a PDE with the Caputo time-fractional derivative of order α, 0 < α ≤ 1 and the fractional spatial derivative of order 2α. Thus the quotient of the derivatives orders is one half, i.e., exactly the same as for the conventional diffusion equation. The case of the two-dimensional α-fractional diffusion equation was considered in Reference [17]. In the n-dimensional case, the entropy of the processes governed by the α-fractional diffusion equation was discussed in Reference [18]. The rest of the paper is organized as follows. In Section 2, we first remind the readers of derivation of the fundamental solution to the n-dimensional time-space-fractional PDE with the Caputo time-fractional derivative of order β, 0 < β < 2 and the fractional spatial derivative (fractional Laplacian) of order α, 0 < α ≤ 2 and then mention its properties that are used in further discussions. Section 3 contains the main results and presents a derivation of the closed form formulas for the Shannon entropy and the entropy production rates of the fractional diffusion processes that are governed by those n-dimensional time-space-fractional PDEs that possess non-negative fundamental solutions (0 < β ≤ 1, 0 < α ≤ 2 for n ∈ N or 1 < β ≤ α ≤ 2 for n = 1). The last section is dedicated to the conclusions.

Fundamental Solution to the Time-Space-Fractional PDE
In this section, for the sake of the reader's convenience we provide a sketch of derivation of the Mellin-Barnes integral representation of the fundamental solution to the time-space-fractional PDE with the Caputo time-fractional derivative and the fractional Laplacian (see Reference [19] for details). This representation is valid for the derivatives orders β ∈ (0, 2), α ∈ (0, 2] [20]. However, our derivation method works only in the case β ∈ (0, 2), α ∈ (1, 2] and thus we first suppose these conditions to hold true and consider the multi-dimensional time-space-fractional PDE in the following form: where D β t is the Caputo time-fractional derivative of the order β and (−∆) α 2 is the fractional Laplacian. The Caputo time-fractional derivative of order β > 0 is defined by the formula where I γ t is the Riemann-Liouville fractional integral given by For a sufficiently well-behaved function f : R n → R and for 0 < α < m, m ∈ N and x ∈ R n , the fractional Laplacian can be represented as a hypersingular integral [21]: with the suitably defined finite differences operator ∆ m h f (x) and the normalization constant d n,m (α). The operator ∆ m h f (x) can be chosen either in the non-centered form or in the centered form: The normalization constant d n,m (α) is given by the formula in the case of the non-centered difference operator and in the case of the centered difference operator. The representation (3) of the fractional Laplacian does not depend on m, m ∈ N provided α < m and is valid with the centered differences operator for all α > 0 and even m or with the non-centered differences operator for all α > 0 except of the case α = 1, 3, 5 . . . , 2[m/2] − 1.
It is worth mentioning that the fractional Laplacian (−∆) α 2 can be also interpreted as a pseudo-differential operator with the symbol |κ| α ( [21,22]): where (F f )(κ) is the Fourier transform of a function f at the point κ ∈ R n defined by In what follows we consider an initial-value problem for the Equation (1) with the Dirichlet initial conditions. For 0 < β ≤ 1, we pose an initial condition in the form If 1 < β < 2, the second initial condition is added: Because the initial-value problem (1), (6) (or (1), (6)- (7), respectively) is a linear one, its solution can be represented in the form where G α,β,n is the fundamental solution, i.e., the solution to the problem (1), (6) with the initial condition or to the problem (1), (6)-(7) with the initial conditions for 0 < β ≤ 1 or 1 < β < 2, respectively, with δ being the Dirac delta function.
To derive a close form formula for the fundamental solution, we apply the Fourier transform (5) with respect to the spatial variable to the Equation (1) and to the initial conditions (6) (7) (in the case β > 1) and get the fractional ordinary differential equation (ODE) and the initial conditionsĜ α,β,n (κ, 0) = 1 (9) in the case 0 < β ≤ 1 or the initial conditionŝ In both cases, the unique solution to the initial-value problem for the fractional ODE (8) with the initial conditions (9) or (10), respectively, has the form [23]: where the Mittag-Leffler function E β (z) is defined as follows: Because of the asymptotic formula [24] the right-hand side of the Equation (11) is from L 1 (R n ) as soon as α > 1. Thus for 1 < α ≤ 2, we can apply the inverse Fourier transform to the Equation (11) and get the representation Because the function E β −|κ| α t β is a radial function, we can rewrite the Equation (14) in the form [21] G α,β,n (x, t) = |x| 1− n 2 (2π) whenever the integral at the right-hand side of (15) converges absolutely or at least conditionally. In this formula, J ν denotes the Bessel function with the index ν.
In the case |x| = 0 (x = (0, . . . , 0)), the formula (14) takes the form that can be rewritten as due to the well-known formula for the multi-dimensional integrals of the radial functions [21]. The integral at the right-hand side of (16) is convergent under the condition 0 < n < α. Thus the fundamental solution G α,β,n is finite at the point |x| = 0 only in the one-dimensional case (remember the condition 1 < α ≤ 2) and has an integrable singularity for other dimensions.
Let us now consider the case x = 0 and apply the variables substitution u = τ α t β in the integral at the right-hand side of the integral representation (15) to get the expression with the auxiliary function This representation will be used in the further discussions of the properties of G α,β,n . Now we employ the technique of the Mellin integral transform to deduce a Mellin-Barnes representation of the fundamental solution G α,β,n (x, t). Recall that the Mellin integral transform of a function f is defined by the formula and the inverse Mellin integral transform has the form If we denote by M ←→ the juxtaposition of a function f with its Mellin transform f * then the convolution theorem for the Mellin integral transform reads as As we can see, the integral at the right-hand side of the Equation (15) can be interpreted as the Mellin convolution of the functions the Mellin integral transform of the Bessel function [25] J ν (2 √ τ) , − (ν/2) < (s) < 3/4, and some elementary properties of the Mellin integral transform [25,26], we arrive at the formulas: Then the Mellin convolution theorem (21) and the Equation (20) for the inverse Mellin transform provide us with the following Mellin-Barnes integral representation of the fundamental solution G α,β,n : where n 2 − 1 2 < γ < min(α, n). The variables substitution s → n − s in the integral at the right-hand side of (22) leads to an equivalent representation with max(n − α, 0) < γ < n.
Comparing the last formula with the expression (17), we get the following Mellin-Barnes representation of the auxiliary function L α,β,n : where the parameter γ satisfies the inequalities max(n − α, 0) < γ < n. Thus we can determine the Mellin integral transform of the function L α,β,n : It is worth mentioning that both the representation (17) and the Equation (25) will play a decisive role for derivation of the closed form formulas for the Shannon entropy and the entropy production rate of the fractional diffusion processes. In particular, the representation (17) shows that the fundamental solution G α,β,n depends on the similarity variable z = |x| 2t β α that in fact determines the form of the entropy production rate of the fractional diffusion processes.
Another important point is the non-negativity property of the fundamental solution. The time-space-fractional PDE governs a fractional diffusion process if and only if its fundamental solution is non-negative and can be interpreted as a spatial probability density function (pdf) evolving in time. As recently shown in [6], these conditions are satisfied for an arbitrary dimension n ∈ N if 0 < β ≤ 1, 0 < α ≤ 2 and additionally for 1 < β ≤ α ≤ 2 in the one-dimensional case. In the next section, we restrict ourselves only to these cases and derive the explicit formulas for the Shannon entropy and for the entropy production rate of the fractional diffusion processes governed by the corresponding time-space-fractional PDEs.
For further properties and numerous particular cases of the fundamental solution G α,β,n we refer the interested readers to [19].

The Entropy Production Rates of the Fractional Diffusion Processes
In this section, we treat the fundamental solution G α,β,n as a probability density function (with respect to the spatial variables) evolving in time (as mentioned above, it is the case for an arbitrary dimension n ∈ N if 0 < β ≤ 1, 0 < α ≤ 2 and additionally for 1 < β ≤ α ≤ 2 in the one-dimensional case) and in particular derive the closed form formulas for the entropy and the entropy production rate of the underlying stochastic processes.
For an n-dimensional continuous random variable evolving in time with the probability density function p(x, t), x ∈ R n , t > 0, the Shannon entropy is defined by the formula and the entropy production rate by its derivative The entropy production rate is an important characteristic of a stochastic process that provides a measure of its irreversibility. For the conventional diffusion process (fundamental solution to the n-dimensional diffusion equation) the Shannon entropy is given by the formula The entropy production rate of the n-dimensional diffusion process is obtained by differentiation of the right-hand side of the Equation (29): The entropy production rate R(G 2,1,n , t) is strictly positive for any t > 0 and thus the corresponding diffusion process can be classified as irreversible. Moreover, R(G 2,1,n , t) decreases with time and goes to zero for t → +∞ that can be interpreted in the sense that the diffusion process tends to a uniform distribution of the particles for t → +∞. Now let us proceed with calculation of the Shannon entropy of the fractional diffusion processes governed by the multi-dimensional time-space-fractional diffusion Equation (1) under the restrictions on the orders of the fractional derivatives that ensure the non-negativity of its fundamental solution (see the previous section). To do this, we employ the representation (17) of its fundamental solution in terms of the auxiliary function L α,β,n .
Substituting (17) into (26) and using the well-known formula for the multi-dimensional integrals of the radial functions ( [21]) we obtain the following chain of equalities: In the last two integrals, we employ the variables substitution τ 2t β/α → τ and arrive at the formula that can be represented in the form where the (time-independent) constants A α,β,n and B α,β,n are given by the expressions in terms of the auxiliary function L α,β,n . The representation (32) allows us to immediately calculate the entropy production rate of the n-dimensional fractional diffusion process governed by the time-space-fractional diffusion equation (1): where the constant A α,β,n is given by (33). Surprisingly, the integral at the right-hand side of the Equation (33) can be calculated in explicit form. Indeed, this integral can be interpreted as the Mellin integral transform of the auxiliary function L α,β,n at the point s = n. However, the Equation (25) for the Mellin integral transform of L α,β,n was derived under the condition n − min(α, n) < Res < n and cannot be directly used at the point s = n. Thus we calculate the integral at the right-hand side of the Equation (33) as the limit of the right-hand side of the Equation (25) as s → n and get the following chain of equalities: The second to last equality is a simple consequence of the formula for the residual of the Gamma-function at the point s = 0: res s=0 Γ(s) = 1. , we arrive at the following formula for the constant A α,β,n : Thus the entropy production rate of the n-dimensional fractional diffusion process governed by the time-space-fractional diffusion Equation (1) can be represented in the following simple form: As in the case of the conventional diffusion, the entropy production rate of the n-dimensional fractional diffusion process is non-negative and decreasing with time. Moreover, in the general case (for any dimension n ∈ N), the entropy production rate is an increasing function of the order of the time-fractional derivative β for β ∈ (0, 1] as expected from the physical and probabilistic interpretations of the slow and conventional diffusion. Because β is restricted to the interval (0, 1] (with the only exception for n = 1), there exists no entropy production rate paradox.
Another important finding follows from comparison of the entropy production rates of the n-dimensional conventional diffusion (Equation (30)) and of the n-dimensional fractional diffusion process (Equation (37)). They are equal in the case i.e., in the case of the α-fractional diffusion equation that was analyzed in detail in [16][17][18]. Thus the α-fractional diffusion equation can be treated as a "right fractionalization" of the conventional diffusion equation.
As one can see from the Equation (37), the entropy production rates and thus other properties of the fractional diffusion processes depend on the quotient β α of the orders of the time-and space-fractional derivatives. Of course, it is a direct consequence of the form of the similarity variable in the representation (17) of the fundamental solution G α,β,n in terms of the auxiliary function L α,β,n . From the other hand, the expression for the similarity variable is determined by the symmetries of the time-space-fractional PDE (1) with respect to the group of its scaling transformations (see [27,28] for discussion of the relevant particular cases and [29] for the general theory of the group invariants for the fractional PDEs).

Conclusions
In this paper, we derived and analyzed the fundamental solution to the n-dimensional time-space-fractional PDE with the Caputo time-fractional derivative of the order β, 0 < β < 2 and the fractional spatial derivative (fractional Laplacian) of the order α, 0 < α ≤ 2. This equation governs a fractional diffusion process if and only if its fundamental solution is non-negative and can be interpreted as a spatial pdf evolving in time that is the case for an arbitrary dimension n ∈ N if 0 < β ≤ 1, 0 < α ≤ 2 and additionally for 1 < β ≤ α ≤ 2 if n = 1. In all these cases, we derived an explicit formula for the entropy production rate of the fractional diffusion processes. It turned out that the rate depends on the orders β and α of the fractional time and spatial derivatives and on the dimension n and is given by the expression β n α t , t being the time variable. The entropy production rate is an increasing function in β. However, one cannot speak about any entropy production paradoxes related to the fractional diffusion processes because the time-space-fractional PDE governs a fractional diffusion process in all dimensions only under the condition 0 < β ≤ 1, i.e., only in the cases of the slow and the conventional diffusion.
Another important finding is that the entropy production rate of the α-fractional diffusion equation (fractional PDE with the Caputo time-fractional derivative of order α, 0 < α ≤ 1 and the fractional spatial derivative of order 2α) is equal to one of the conventional diffusion equations. Thus the α-fractional diffusion equation can be treated as a "right fractionalization" of the conventional diffusion equation.
It would be important to study the problems discussed in this paper for other types of the fractional PDEs. In particular, the fractional diffusion equations with the Riesz-Feller fractional spatial derivative and/or the Riemann-Liouville time-fractional derivative would be worth investigating.