{'fraction_non_alphanumeric': 0.055419934555653666, 'fraction_numerical': 0.013348569054173375, 'mean_word_length': 4.155576382380506, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Some recent advances in theory and simulation of fractional diffusion processes

To offer a view into the rapidly developing theory of fractional diffusion processes we describe in some detail three topics of present interest: (i) the well-scaled passage to the limit from continuous time random walk under power law assumptions to space-time fractional diffusion, (ii) the asymptotic universality of the Mittag-Leffler waiting time law in time-fractional processes, (iii) our method of parametric subordination for generating particle trajectories.


Introduction
The field of fractional (more generally anomalous) diffusion processes in recent decades has won more and more interest in applications in the sciences, in physics and chemistry, and even in finance. Here we will give some insight into this rapidly developing field. Instead of trying to cover the whole spectrum of this field we will focus on three contrasting aspects to which we ourselves have affinity by our research, thereby trying to keep this paper as self-contained as possible. We hope not only to whet the appetite among people not fully familiar with the subject but also to propagate our methodological viewpoints in the fractional diffusion community.
Viewing fractional diffusion processes as generalization of the familiar classical diffusion process, we begin by considering the Cauchy problem for the classical diffusion equation In fractional diffusion equations (see Section 2 for explanations) the differentiations with respect to t and x are replaced by differentiations of non-integer order. For equation (1.1) it is well known that its solution u(x, t) has the essential properties we expect from a diffusion process, that is a process of re-distribution in space x and time t. Considering u(x, t) as the spatial density of an extensive quantity, e.g. mass, charge, or probability, we have (a), (b) and (c): (a) conservation of the total quantity: +∞ −∞ u(x, t) dx = +∞ −∞ u(x, 0 + ) dx, ∀ t > 0. (b) preservation of non-negativity: u(x, 0 + ) ≥ 0, ∀ x ∈ R implies u(x, t) ≥ 0, ∀ x ∈ R, t > 0. (c) Another essential characteristic of problem (1.1) concerns the law of spreading (or dispersion) of the quantity. With the special initial condition u(x, 0 + ) = δ(x) (the Dirac generalized function), the variance grows linearly in time, that is σ 2 (t) := diffusion for processes having properties (a) and (b) but not necessarily (c), we follow the custom to call processes of anomalous diffusion those in which, for initial condition u(x, 0 + ) = δ(x), the variance does not exhibit essentially linear grow with t → ∞. Among these processes we single out the sub-diffusive ones for which the variance grows (for large t) more slowly than linearly,and the super-diffusive for which it grows (for large t) faster than linearly, or even does not exist (i.e is infinite). Focusing our attention to the space-time fractional diffusion equation (i.e. a generalization of the classical diffusion equation (1.1) via suitable pseudo-differential operators interpreted as time and space derivatives of fractional order), we will discuss the construction and properties of its fundamental solution (obtained for f (x) = δ(x)) and its approximation by continuous time random walk. This generalized diffusion equation has found wide interest among researchers in recent 20 years. After a general survey of basic facts we will go into details of three distinct but related topics. As this paper cannot be a substitute for an extensive monograph, our presentation will naturally be biased by our and our close co-workers' contributions. We will meet the two complementary aspects of diffusive processes. The first is the macroscopic aspect: the structure of the fundamental solutions in dependence on space x and time t, in particular their scaling properties and asymptotics. Here u(x, t) is viewed as the density with respect to x evolving in t. The second is the microscopic aspect: the structure of the trajectories (paths) of particles subject to such process. Here u(x, t) is viewed as the sojourn probability density with respect to x evolving in t.
The rough structure of our paper is as follows. In Section 2 we will give a survey of the space-time fractional diffusion equation and the essential properties of its fundamental solution. Sections 3 and 4 are devoted to topic (i): continuous time random walk (CTRW), Section 5 to topic (ii): power laws and well-scaled passage to the diffusion limit, Section 6 to topic (iii): our method of parametric subordination for exact simulation of trajectories. Finally, in Section 7, we will draw some conclusions.
Of course, there are more significant advances than we can report here. Let us only mention processes with distributed orders of fractional differentiation and multi-dimensional processes and cite, e.g. the papers [1,5,6,46] and the fundamental monograph by Meerschaert & Scheffler [33]. Much work has been done in numerical analysis of difference schemes for calculating densities of fractional diffusion processes, see e.g. [22]. We apologize to all authors who feel that we appreciate their work insufficiently. Throughout this paper we will make liberal use of generalized functions in the sense of Gel'fand and Shilov [10], so avoiding the cumbersome notations of measures and functionals, and usually we will assume the occurring (generalized) functions so well-behaved that our manipulations are allowed. By this we hope our presentation to be welcome to researchers in applications as well as inspiring for pure mathematicians looking around outside the ivory tower.

The space-time fractional diffusion
We begin by considering the Cauchy problem for the (spatially onedimensional) space-time fractional diffusion equation where α, θ, β are real parameters restricted to the ranges Here t D β * denotes the Caputo-Dzherbashyan fractional derivative of order β, acting on a sufficiently well-behaved function f (t) of the time variable t, and x D α θ denotes the Riesz-Feller fractional derivative of order α and skewness θ, acting on a sufficiently well-behaved function g(x) of the space variable x, In the symmetric case θ = 0, which later will be of our main interest, formula (2.4a) simplifies to (2.4 ′ a) The above representations of the space fractional derivatives are based on a suitable regularization of hyper-singular integrals.
In the limits β = 1 and α = 2 (so θ = 0) we recover the first time derivative df (t) dt and the second space derivative d 2 g(x) dx 2 , respectively. These representations mirror the fact that time-fractional (for 0 < β < 1) processes are processes with long memory (see for this also Section 5) whereas space fractional (for 0 < α < 2) are processes with spatial longrange interactions. For more information on these operators we refer to [14,15,26,28,38,40]. Let us note that the solution u(x, t) of the Cauchy problem (2.1), known as its Green function or fundamental solution, can be viewed as the density of an extensive quantity or as the probability density in the spatial variable x, evolving in time t. In the case α = 2 (hence θ = 0) and β = 1 we recover the standard diffusion equation for which the fundamental solution is the Gaussian density with variance σ 2 = 2t.
Writing, with Re[s] > σ 0 , κ ∈ R, the transforms of Laplace and Fourier as We will freely use the convolution theorems pertinent to these transforms, defining for generic functions: x ∈ R, and the convolution powers f * n (t) and g * n (x) as n-fold convolutions ( n ≥ 0). Note that f * 0 (t) = δ(t), g * 0 (x) = δ(x). For mathematical details we cite [14,38] on the Caputo-Dzherbashyan derivative and [40] on the Feller potentials. For the general theory of pseudo-differential operators and related Markov processes the interested reader is referred to the excellent volumes by Jacob [23]. Let us here recall the representation in the Fourier-Laplace domain of the (fundamental) solution of (2.1). Using δ(κ) ≡ 1 we have from (2.1) s β u(κ, s) − s β − 1 = −|κ| α i θ sign κ u(κ, s) , For explicit expressions and plots of the fundamental solution of (2.1) in the space-time domain we refer the reader to [28]. There, starting from the fact that the Fourier transform u(κ, t) can be written as a Mittag-Leffler function with complex argument, the authors have derived a Mellin-Barnes integral representation of u(x, t) with which they have proved the non-negativity of the solution for values of the parameters {α, θ, β} in the range (2.2) and analyzed the evolution in time of its moments. In particular for {0 < α < 2, β = 1} we obtain the densities of the stable processes of order α and skewness θ. The representation of u(x, t) in terms of Fox H-functions can be found in [29]. We note, however, that the solution of the space-time fractional diffusion equation (2.1) and its variants has been investigated by several authors as pointed out in the bibliography in [28]: here we refer to some of them, [2,3,31,32,34], where the connection with the CTRW was also pointed out. Henceforth our attention, if not said explicitly otherwise, will be focussed on the symmetric case θ = 0. In this case u(x, t) is an even function of x and we get from (2.6) the Fourier-Laplace representation which allows us demonstration of the conservation property (a), namely +∞ −∞ u(x, t) dx ≡ 1 for all t > 0, and calculation of the variance σ 2 (t) = x 2 (t) (the quadratic measure of the spread as function of t) of a diffusing particle. From (2.7), more generally already from (2.6), we find by aid of well-known properties of the Fourier transform u(0, s) = 1/s, hence For the variance we find . (2.8) Using the Mittag-Leffler function see e.g. [14,38], we find by Laplace inversion the convergent series from which for t > 0 we get (2.10) Admitting that these calculations are formal we leave the task of strict justification to pure mathematicians.

The continuous time random walk
In the Sixties and Seventies of the past century Montroll and Weiss and Scher (just to cite these pioneers) published a series of papers for modelling rather general diffusion processes by random walks on lattices, see e.g. [35,36], and the book by Weiss [48] with references therein. Initiated by their activity the concept of continuous time random walk became popular in physics. CTRWs are good phenomenological models for several types of diffusion, in particular the microscopic aspects of a particle jumping from point to point admitting waiting times between jumps. Allowing all space instead of only a lattice, a CTRW can be viewed as a compound renewal process or a renewal process with reward (see [7]), or a random walk subordinated to a renewal process. Let us recall the basic notions of the CTRW theory.
A CTRW is generated by a sequence of independent identically distributed (iid) positive random waiting times T 1 , T 2 , T 3 , . . . , each having the same probability density function φ(t) , t > 0 , and a sequence of iid random jumps X 1 , X 2 , X 3 , . . . , in R , each having the same probability density w(x) , x ∈ R . Setting t 0 = 0 , t n = T 1 + T 2 + . . . T n for n ∈ N , the wandering particle makes a jump of length X n in instant t n , so that its position is x 0 = 0 for 0 ≤ t < T 1 = t 1 , and x n = X 1 + X 2 + . . . X n , for t n ≤ t < t n+1 .
We require the distribution of the waiting times and that of the jumps to be independent of each other. We allow the probability densities φ and w to be generalized functions in the sense of Gel'fand and Shilov [10].
Natural probabilistic arguments lead us to the integral equation for the probability density p(x, t) (a density with respect to the variable x) of the particle being in point x at instant t , see e.g. [19,30,42,43,44], dt ′ denotes the probability that at instant t the particle is still sitting in its starting position x = 0 . Clearly, (3.1) satisfies the initial condition p(x, 0) = δ(x), and p in place of u has the properties (a) and (b) of Section 1. In the Fourier-Laplace domain Eq. (3.1) reads as and using Ψ(s) = (1 − φ(s))/s , explicitly This representation is known as the the Montroll-Weiss equation, so named after the authors of [36], who derived it in 1965 as the basic equation for the CTRW. By inverting the transforms one can find the evolution of the sojourn density p(x, t) for time t running from zero to ∞. In fact, recalling that | w(κ)| < 1 and | φ(s)| < 1, if κ = 0 and s = 0, Eq. (3.3) becomes and we promptly obtain the series representation of the continuous time random walk, see e.g. [7] (Ch. 8, Eq. (4)) or [48] (Eq.(2.101)), where the functions v n (t) and w n (x) are obtained by iterated convolutions in time and in space, v n (t) = (Ψ * φ * n )(t), and w n ( In the R.H.S of Eq (3.5) we have isolated the first singular term related to the initial condition p(x, 0) = Ψ(0) δ(x) = δ(x). The representation (3.5) can be found without detour over (3.1) by direct probabilistic reasoning and transparently exhibits the CTRW as a subordination of a random walk to a renewal process: it can be used as 8 starting point to derive the Montroll-Weiss equation, as it was originally recognized by Montroll and Weiss [36].
With the special choice φ(t) = me −mt , m > 0, equation (3.1) describes the compound Poisson process. It reduces after some manipulations (best carried out in the Fourier-Laplace domain) to the Kolmogorov-Feller equation: Then from (3.5) we obtain the series representation Note that only in this case the corresponding stochastic process is Markovian.

Relevance of power laws and well-scaled passage to the diffusion limit
In this section we work out the effect the power laws for the distributions of jumps and waiting times have on the limiting properties of continuous time random walks. By power law we usually mean a law of decaying at infinity like a power (with negative exponent) of the relevant independent variable. There is a vast amount of literature on such laws. Let us recommend just two items, namely [37] and [45]. What we are going to show now is that appropriate power laws for jumps and waiting times are microscopic models for fractional diffusion processes.
To be sufficiently general we introduce the cumulative functions for waiting times and jumps. With our densities φ(t) and w(x) we define These functions are non-decreasing, As they may have points of discontinuity we agree on the provision that equations in which they occur are meant to hold at points of continuity. The notion of power law concerns the behaviour of Φ(t) and W (x) for large t and large |x|, roughly the mode of decrease of the tails Not wanting to overload our presentation we assume spatial symmetry (θ = 0) with respect to x = 0 and avoid decoration of asymptotic behaviours with slowly varying functions. For such neglected generalizations we refer to [13,17]. The parameters α and β of equation (2.1) play essential roles yielding power laws in the strict sense if 0 < α < 2, 0 < β < 1, but still formal analogies in the extreme cases α = 2, β = 1.
The question of interest is the long-time, wide-space behaviour of a CTRW under power law assumptions for waiting times and jumps, i.e. the appearance of such CTRW when observed after long time and from far away. To bring the distant future and the far-away space into near sight we multiply time intervals by a small positive factor τ , space intervals by a small positive factor h, so making large intervals numerically of moderate size, intervals of moderate size numerically small. Essentially this means changing the units of time and space from 1 to 1/τ and 1/h, respectively. We then obtain the asymptotic behaviour by sending τ and h to zero, in a specially combined way, namely under the requirement of a scaling relation, honouring which we call the whole procedure "well-scaled passage to the diffusion limit". In this limit we will obtain a process obeying the spacetime fractional diffusion equation (2.1) with θ = 0 , a fact that justifies the CTRW approach. Conversely, we can consider a CTRW as a model "in the small" of a space-time fractional diffusion process. For generalization to skewed processes (θ = 0) see [17].
As we carry out the essential work in the Fourier-Laplace domain we use the fact that the behaviour of functions f (t), g(x) in the infinite is mirrored in that of their transforms f (s), g(κ) for s → 0 + , κ → 0. The lemmata we need are provided by the Tauber-Karamata theory for which we refer to [4] and [8], they can also be distilled from the Gnedenko theorem on the domains of attraction of stable probability laws, see [11]. See [13] and [17] for more general versions. What we need is contained in the following two lemmata which are simplified versions of more general facts. For the reader's convenience we give proofs in the Appendix.
Then we have the asymptotics Then we have the asymptotics Assuming the conditions of these two lemmata satisfied we are ready for passing to the diffusion limit. We multiply the jumps X k by a factor h > 0, the waiting times T k by a factor τ > 0. So we get a transformed random walk S n (h, τ ) = n k=1 hX k with jump instants t n (τ ) = n k=1 τ T k that we now investigate with the aim of passing to the limit h → 0, τ → 0 under a scaling relation between h and τ yet to be established. As it is convenient to work in the Fourier-Laplace domain we note that the density φ τ (t) of the reduced waiting times τ T k and the density w h (x) of the reduced jumps The corresponding transforms are simply φ τ (s) = φ(sτ ), w h (κ) = w(κh). We are interested in the sojourn probability density p h,τ (x, t) of the particle subject to the transformed random walk. In analogy to the Montroll-Weiss equation Considering now s and κ fixed and = 0 we find for h → 0, τ → 0 from the Master Lemmata (replacing there κ by κh, s by sτ ) by a trivial calculation, omitting asymptotically negligible terms, the asymptotics (4.3) with (4.4).
So we see that, for every fixed real κ = 0 and positive s, as h and τ tend to zero under the scaling relation r(h, τ ) ≡ 1. Comparing with (2.7) we recognize here u(κ, s) as the combined Fourier-Laplace transform of the solution to the Cauchy problem (2.1) with θ = 0. Invoking now the continuity theorems of probability theory (compare [8]), bypassing some analytical subtleties, we see that the time-parametrized sojourn probability density converges weakly (or in law) to the solution of the Cauchy problem (2.1) with θ = 0. This weak convergence can be taken as justification for approximate simulation of trajectories (particle paths) by CTRW's with power law jumps and waiting times so chosen that routines are available for producing the needed random numbers (see e.g. [13]).

The Mittag-Leffler waiting time law and timefractional processes
We now offer another view to the transition from CTRW (under power law assumptions) to fractional diffusion, separating the temporal and spatial limiting procedures, thereby getting among other matters the time-fractional CTRW discussed in the pioneering paper [21].
Turning our attention to time-fractional processes we present in a condensed way some results from our papers [18,19,30]. We will see the importance of the Mittag-Leffler waiting time density φ M L (t) and the corresponding survival function Ψ M L (t) with their Laplace transforms displayed here: with E β defined in (2.9). Throughout we assume 0 < β ≤ 1. We recall that for β = 1 we recover φ M L (t) = exp(−t), Ψ M L (t) = exp(−t).
The importance of these functions cannot be overestimated, they also play an essential role in the theory of fractional relaxation (see e.g. [14]). The relevance of the Mittag-Leffler waiting time law for time-fractional CTRW has been put in bright light by Hilfer and Anton in 1995 [21]. Fulger, Scalas and Germano [9] pay special attention to its use as waiting time law in CTRW simulation . In the Sixties of the past century it has been found by Gnedenko and Kovalenko [12] as a limiting law in thinning (rarefaction) of a renewal process under the power law assumptions of our Master Lemma 2, but they only gave its Laplace transform without identifying it as a function of Mittag-Leffler type. We will show, without invoking the concept of thinning, that it is a universal limiting law for long-time behaviour of a renewal process under a power law regime.
For a general CTRW we will show how via the concept of the memory function we can separate the passages of the scaling factors τ and h (of the preceding Section 4) to zero, thus avoiding the simultaneous use of the continuity theorems for the transforms of Laplace and Fourier.

Manipulations: rescaling and respeeding
To introduce the memory function H(t) and explain its meaning we recall the CTRW theory of Section 3, in particular equations (3.1) and (3.3). We need this general theory already for embedding into it the renewal theory.
In fact: we can view a pure renewal process as a special type of CTRW, namely one in which all jumps are of fixed size 1 and the position x of the wandering particle in space plays the role of the counting number n of the renewal process. However, we prefer to continue working in the general CTRW context. Introducing formally in the Laplace domain the auxiliary function and assuming that its Laplace inverse H(t) exists, we get, following [30], in the Fourier-Laplace domain the equation and in the space-time domain the generalized Kolmogorov-Feller equation with p(x, 0) = δ(x). With the special choice of the power law memory function whose Laplace transform is For 0 < β < 1 we have the time-fractional Kolmogorov-Feller equation Let us now consider two types of manipulations on the CTRW: A: rescaling, B: respeeding.
(A) means, as in Section 4, change of the unit of time. With the positive scaling factor τ we replace the waiting time T by τ T , again intending τ << 1. In a moderate span of time we will so have a large number of jump events. Again we get the rescaled waiting time density and its corresponding Laplace transform as φ τ (t) = φ(t/τ )/τ , φ τ (s) = φ(τ s). By decorating also the density p with an index τ we obtain this rescaled CTRW integral equation in the Fourier-Laplace domain as fixing s as required by the continuity theorem of probability for Laplace transforms, the asymptotics φ(s) ∼ 1 − λs β for 0 < τ → 0 implies corresponding to the density φ M L (t). This formula expresses the asymptotic universality of the Mittag-Leffler waiting time law that includes the exponential law for β = 1. It says that our general power law waiting time density is gradually deformed into the Mittag-Leffler waiting time density.
under combined rescaling with τ and respeeding with a = λτ β . This attraction property of the Mittag-Leffler waiting time distribution with respect to power law waiting times (with 0 < β < 1) is a kind of analogy to the attraction of sums of power law jump distributions by stable distributions.

Diffusion limit in space
We can obtain from (5.5) the fractional Kolmogorov-Feller equation (5.10) for time fractional CTRW by direct insertion of the Mittag-Leffler memory function into the equation or, as the previous considerations show, by manipulating it under power law assumption for the waiting time and passing to the limit. We have not yet operated on the jumps. To do this now, we assume the conditions of Master Lemma 1 to hold. Then, by another respeeding, in fact an acceleration (that we earlier had carried out in [19]), we will arrive at diffusion processes fractional in space. We have three choices: Note that (B) is just a special case of (A) but of particular relevance.In all three cases we rescale the jump density by a factor h > 0, replacing the random jumps X by hX. This means changing the unit of measurement in space from 1 to 1/h, with 0 < h ≪ 1, so bringing into near sight the far-away space. The rescaled jump density is w h (x) = w(x/h)/h, corresponding to w h (κ) = w ( hκ). still with H(s) as in (5.3). Here −|κ| α is the symbol of the Riesz pseudodifferential operator x D α 0 (known as the Riesz fractional derivative of order α) obtained from the Riesz-Feller fractional derivative for θ = 0, see (2.1) and (2.4). We thus arrive at the integro-pseudo-differential equation (5.25) Comments: By this rescaling and acceleration the jumps become smaller and smaller, their number in a given span of time larger and larger, the waiting times between jumps smaller and smaller. In the limit there are no waiting times anymore, the original waiting time density φ(t) is now only spiritual, but still determines via H(t) the memory of the process. Eq. Comments: Let us point out an advantage of splitting the passages τ → 0 and h → 0. Whereas by the combined passage as in choice (C), if done in the well-scaled way (5.26), the mystical concept of respeeding can be avoided, there arises the question of correct use of the continuity theorems of probability. There is one continuity theorem for the Laplace transform, one for the Fourier transform, see [8]. Possible doubts whether their simultaneous use is legitimate vanish by applying them in succession, as in our two splitting methods. For a more detailed discussion of mathematical aspects and the involved stochastic processes we refer to our recent paper [18].

Subordination in stochastic processes
The common method for simulating particle trajectories consists in interpreting the concept of subordination (see [8]) as one of transforming a stochastic process Y (t * ) (we call it the parent process) where t * is not the physical but an operational time into a process X(t) in physical time t, by generating the operational time t * from the physical time t via a positively oriented stochastic process T * (t) to arrive at the representation X(t) = Y (T * (t)). For simulation one then needs a routine for generating the process T * (t). For simulating trajectories in space-time fractional diffusion (2.1) this requires simulation of the hitting time process which is inverse to the stable subordinator in Feller's parametrization [8], the stable process of order β and skewness −β.
There are routines available for simulating stable variates, see e.g. [24,25]. But we do not know of easy routines for inverting a stable subordinator. Therefore, we have developed our method of parametric subordination which, by starting from the operational time t * allows construction of trajectories by producing: FIRST an α-stable process x = Y (t * ) with skewness θ for for the position x in space, SECOND an extreme positive-oriented stable process t = T (t * ) of order β with skewness −β that we call the leading process. Then we get, in the (t, x)-plane, the parametrized graph t = T (t * ), x = Y (t * ) of a desired trajectory of the process X(t) corresponding to eq. (2.1) as For the more general situation we refer to our recent paper [20]. This method is exact in the sense that it allows us to produce snapshots of a true particle path. Let us in this context also draw attention to the recent paper [27] by Kleinhaus and Friedrich.
Let us sketch how this method directly arises from the CTRW approximation under appropriate power law assumptions on the waiting time and the jumps. Compare also with [39] and [47]. For mathematical details we refer to [20]. We start with the equations (3.4) and (3.5), and rescale waiting times and jumps again with positive factors τ and h. In the Fourier-Laplace domain we assume, a bit more general than in our Master lemmata, the power-law conditions We see that we have a combination of two Markov processes happening in discrete time, one giving a jump (with density w(x)) in space x at every instant n where n is the running index, the other one giving a positive jump in time t at every instant n.
Using the effect of the rescaling on Eq. (3.5) and correspondingly decorating it with additional indices we get in the Fourier-Laplace domain Separately we treat the powers φ(τ s) n and ( w(hκ)) n , so avoiding the problematic simultaneous inversion of the diffusion limit from the Fourier-Laplace domain into the physical domain. Observing from (6.1) φ(τ s) n ∼ 1 − λ(τ s) β n , we relate the running index n to the presumed operational time t * by n ∼ t * /(λ τ β ), and for fixed s (as required by the continuity theorem of probability theory), by sending τ → 0 we get φ(τ s) n ∼ 1 − λ τ β s β t * /(λτ β ) → exp −t * s β . Here the Laplace variable s corresponds to physical time t, and in Laplace inversion we must treat t * as a parameter. Hence, in physical time exp(−t * s β ) corresponds tō with ḡ β (s) = exp(−s β ). Hereḡ β (t, t * ) is the totally positively skewed stable density (with respect to the variable t) evolving in operational time t * according to the "space"-fractional equation where t plays the role of the spatial variable. Analogously, observing from (6.2) ( w(hκ)) n ∼ 1 − µ(h|κ|) α i θsign κ n , with the aim of obtaining a meaningful limit we set n ∼ t * /(µh α ), and find, by sending h → 0 + , the relation the Fourier transform of a θ-skewed α-stable density f α,θ (x, t * ) evolving in operational time t * . This density is the solution of the space-fractional equation The two relations between the running index n and the presumed operational time t * require the (asymptotic) scaling relation λ τ β ∼ µ h α , that for purpose of computation we simplify to Replacing t * by t * ,n = nλτ β , using the asymptotic results obtained for the powers φ(τ s) n and ( w(hκ)) n , furthermore noting (1 − φ(τ s))/s ∼ s β−1 λ τ β , we finally obtain from (6.3) the Riemann sum (with increment λτ β ) exp −nλτ β s β + |κ| α i θsign κ λ τ β , (6.8) and hence the integral For the limiting process u β (x, t) this means Observe that the RHS of this equation is just another way of writing the RHS of equation (2.6) which is the Fourier-Laplace solution of the spacetime fractional diffusion equation (2.1). By inverting the transforms we get after some manipulations (compare [31]) in physical space-time the integral formula of subordination Eq. (6.11) is basic for the conventional concept of subordination where there are also two processes involved. One is the unidirectional motion along the t * axis representing the operational time. This motion happens in physical time t and the pdf for the operational time having value t * is (as density in t * , evolving in physical time t) given by (6.12). In fact, by substituting The operational time t * stands in analogy to the counting index n in Eqs. (3.5) and (6.3). The other process is the process described by Eq. (6.6), a spatial probability density for sojourn of the particle in point x evolving in operational time t * , namelyū β (x, t * ) = f α,θ (x, t * ). We get the solution to the Cauchy problem (2.1), namely the pdf u(x, t) = u β (x, t) for sojourn in point x, evolving in physical time t, by averagingū β (x, t * ) with the weight function g β (t * , t) over the interval 0 < t * < ∞ according to (6.11).

Trajectories for space-time fractional diffusion
In the series representation (3.5) of the CTRW the running index n (the number of jumps having occurred up to physical time t) is a discrete operational time, proceeding in unit steps. To this index n corresponds the physical time t = t n , the sum of the first n waiting times, and in physical space the position x = x n , the sum of the first n jumps, see Section 3. We have two discrete Markov processes (discrete in operational time n), namely a random walk in the space variable x, with jumps X n , and another random walk (only in positive direction) of the physical time t, making a forward jump T n at every instant n. Rescaling space and physical time by factors h and τ , observing the scaling relation (6.7), we introduce, by sending h → 0 and τ → 0, the continuous operational time Then in the diffusion limit the spatial process becomes an α-stable process for the position x =x =x(t * ), whereas the unilateral time process becomes a unilateral (positively directed) β-stable process for the physical time t = t =t(t * ). A trajectory of a diffusing particle in physical coordinates can be produced by combining in the (t, x) plane the two Markovian random evolutions x =x =x(t * ) , t =t =t(t * ) , (6.15) obeying the stochastic differential equations (compare with [27]) dx = d(Lévy noise of order α and skewness θ) , dt = d(one sided positively oriented Lévy noise of order β) . (6.16) This gives us in the (t, x) plane the t * -parametrized particle trajectory that by elimination of t * we get as x = x(t). We suggest to call this procedure "construction of a particle trajectory by parametric subordination". Note that the process t = T (t * ) yielding the second random function in (6.16) has the properties of a subordinator in the sense of Definition 21.4 in [41].
6.2 Numerical results for the symmetric case θ = 0 For numerical simulation of trajectories we proceed in three steps. First, let the operational time t * assume N discrete equidistant values in a given interval [0, T ], that is t * ,n = nT /N, n = 0, 1, . . . , N . As a working choice we take T = 1 and N = 10 6 . Then produce N independent identically distributed (iid) random deviates, Y 1 , Y 2 , . . . , Y N having a symmetric stable probability distribution of order α, see the book by Janicki [24] for a useful and efficient method to do that. Now, with the points the couples (t * ,n , x n ), plotted in the (t * , x) plane (operational time, physical space) can be considered as points of a true trajectory {x(t * ) : 0 ≤ t * ≤ T } of a symmetric Lévy motion with order α corresponding to the integer values of operational time t * = t * ,n . In this identification of t * with n we use the fact that our stable laws for waiting times and jumps imply λ = µ = 1 in the asymptotics (6.1) and (6.2) and τ = h = 1 as initial scaling factors in (6.3) and (6.7). In order to complete the trajectory we agree to connect every two successive points (t * ,n , x n ) and (t * ,n+1 , x n+1 ) by a horizontal line from (t * ,n , x n ) to (t * ,n+1 , x n ), and a vertical line from (t * ,n+1 , x n ) to (t * ,n+1 , x n+1 ). Obviously, this is not the 'true' Lévy motion from point (t * ,n , x n ) to point (t * ,n+1 , x n+1 ), but from the theory of CTRW we know this kind of discrete random process to converge in the appropriate sense to Lévy motion. The points (t * ,n , x n ) are points of a true Lévy motion, as can be shown by observing the infinite divisibility and self-similarity of stable laws.
As a second step, we produce N iid random deviates, T 1 , T 2 , . . . , T N having a stable probability distribution with order β and skewness −β (extremal stable distributions). Then, consider the points The final (third) step consists of plotting points (t(t * ,n ), x(t * ,n )) in the (t, x) plane, namely the physical time-space plane, and connecting them as before. So one gets a discrete approximation of a particle trajectory of spatially symmetric (θ = 0) fractional diffusion with parameters α and β.
Now as the successive values of t * ,n and x n are generated by successively adding the relevant standardized stable random deviates, the obtained sets of points in the three coordinate planes: (t * , t), (t * , x), (t, x) can, in view of infinite divisibility and self-similarity of the stable probability distributions, be considered as snapshots of the corresponding true random processes occurring in continuous operational time t * and physical time t, correspondingly. Clearly, fine details between successive points are missing. They are hidden: -In the (t * , x) plane in the horizontal lines from (t * ,n , x n ) to (t * ,n+1 , x n ) and the vertical lines from (t * ,n+1 , x n ) to (t * ,n+1 , x n+1 ).
-In the (t * , t) plane in the horizontal lines from (t * ,n , t n ) to (t * ,n+1 , t n ) and the vertical lines from (t * ,n+1 , t n ) to (t * ,n+1 , t n+1 ). -In the (t, x) plane in the horizontal lines from (t n , x n ) to (t n+1 , x n ) and the vertical lines from (t n+1 , x n ) to (t n+1 , x n+1 ).
The well-scaled passage to the diffusion limit here consists simply in regularly subdividing the {t * } intervals of length 1 into smaller and smaller subintervals (all of equal length τ and adjusting the random increments of t and x according to the requirement of self-similarity, namely taking, respectively, the waiting times and spatial jumps as τ 1/β multiplied by a standard extreme β-stable deviate, τ 1/α multiplied by a standard (in our special case: symmetric) α-stable deviate, respectively, as required by the self-similarity properties of the stable probability distributions). Furthermore if we watch a trajectory in a large interval of operational time t * , the points (t * ,n , x n ) and (t * ,n+1 , x n+1 ) will in the graphs appear very near to each other in operational time t * and aside from missing mutually cancelling jumps up and down (extremely near to each other) we have a good picture of the true processes.

Conclusions
Fractional diffusion processes as models of anomalous diffusion are gaining increasing popularity not only among science researchers but also among more or less pure mathematicians. For the latter they offer fascinating opportunities for applying pseudo-differential operators and other powerful analytic instruments, e.g. those of fractional calculus that in recent decades has made remarkable advances.
In the field of anomalous diffusion we meet challenges for the experimental sciences, for mathematical modelling of real processes and their simulation, for investigation of the underlying evolution processes (the macroscopic aspect) and the fine-structure of their particle trajectories (the microscopic aspect), and for numerical analysis and computational treatment of less common problems. In our presentation we have discussed three topics of current interest that make visible the large arsenal of tools required.
We have completed the proof.

Proof of Master Lemma 2.
In the easy case (A) β = 1 the statement 1 − φ(s) = ρ s is a consequence of the well-known fact ρ = − φ ′ (0). In case (B) 0 < β < 1 we invoke Corollary 8.1.7 of [4]. It says, among other things, that for a probability distribution function Φ(t) vanishing for t < 0 the relation where L(t) is a slowly varying function at infinity, implies the relation 1 − φ(s) ∼ s β L(1/s) for s → 0 + . Now taking L(t) ≡ c Γ(1 − β)/β we get where we have used the reflection formula for the gamma function. The proof is complete. 28