Letters $\mathbf{P}$ and $\mathbf{E}$ will be used everywhere to denote probability and expectation. The indicator function of a set M will be denoted $\mathbf{1}\left(M\right)$.

The CTRWs (continuous time random walks) found numerous applications in physics. The scaling limits of these CTRW were analysed by many authors, see e.g. [1,2,3] and references therein. The crucial points (realised initially by physicists, see [4,5]) were as follows: (i) the limits of scaled CTRWs yield Markov processes time-changed by inverse stable subordinators and (ii) these limiting processes solve fractional in time partial differential equations (PDEs). For general properties of fractional PDEs we can refer to [6,7].

The simplest CTRWs in ${\mathbf{R}}^d$ are specified by an i.i.d. sequence $\left\{W_i\right\}$ of positive random variables (representing waiting times) with a distribution $Q\left(dw\right)$ and an i.i.d. sequence of random variables $\left\{X_i\right\}$ in ${\mathbf{R}}^d$ (representing sizes of jumps) with a distribution $P\left(dx\right)$. This pair of random sequences defines the following process in ${\mathbf{R}}^d$. A particle starting at a point $x\in {\mathbf{R}}^d$ at a time w waits a random time $W_1$ in x, then the particle makes a jump to a random point $x+X_1$, then it waits at $x+X_1$ a random time $W_2$, and then jumps to a random point $x+X_1+X_2$, etc. The total position and waiting time to the time t are therefore respectively, where $\left[t\right]$ denotes the integer part of t. The pair $(X_x\left(t\right),W_w\left(t\right))$ is a Markov chain in ${\mathbf{R}}^d\times \mathbf{R}$ with the transition operator

The CTRW is then defined as the first (spatial) coordinate of this pair subordinated by the inverse process $N_K=sup\{t:W\left(t\right)\le K\}$ to the second coordinate, that is, as the process $x+{\sum }_{j=1}^{N_K}{X}_j$. In physics such processes are usually referred to as Lévy flights. Notice that the process $N_K$ can be equivalently defined by the property that $N_K=n⟺W\left(n\right)\le K<W_{n+1}$ or yet equivalently by the property that $N_K\ge n⟺W\left(n\right)\le K$.

One can naturally extend these simplest CTRWs to the case when the variables $X_i,W_i$ are not independent. For instance, they may have position dependent jumps with dependent $X_i,W_i$. In such cases the transition operator of the corresponding Markov chain takes the form with some transition stochastic kernel Q. Of course, they can also depend on t. The natural general extension leads us to look at the processes of type $X\left(N_K\right)$, where $(X,W)\left(t\right)$ is some Markov process in ${\mathbf{R}}^d\times \mathbf{R}$ with an increasing second coordinate $W\left(t\right)$ and $N_K=sup\{t:W\left(t\right)\le K\}$ is the inverse process to the second coordinate.

Specific new feature of the limiting equations for dependent $(X_i,T_i)$ is that the corresponding fractional derivatives do not separate in time and spatial additive terms, but appear in a certain mixed form. Moreover, as was found in seminal papers [8,9], in case of dependent times and sizes of jumps, it is natural to distinguish two versions of subordination defining CTRWs: $X\left(N_K\right)$ and $X(N_K+1)$, referred to as lagging and leading CTRWs, or, in a different terminology, undershooting CTRWs (or just CTRWs) and overshooting CTRWs.

Popular examples of processes arising from dependent $X_i,T_i$ represent the Lévy walks, which are the main objects of the present study. We refer to [10,11,12,13,14] for recent results and up-to-date reviews of mathematical and physical literature on these processes.

When moving according to a simple Lévy walk, the particle, instead of waiting a random time $W_i$ in a certain location $X_{i-1}$, is moving during the time $W_i$ with some constant velocity $V_i$ chosen from a given distribution, thus arriving at $X_i=X_{i-1}+V_i{W}_i$. In the simplest case one assumes the sequences of pairs $\left\{(V_i,W_i)\right\}$ to be i.i.d. with also independent $V_i$ and $W_i$.

Strictly speaking, Lévy walks are quite different from Lévy flights, as the former have continuous trajectories (particles move at constant velocities), while Lévy flights are jump-type processes with discontinuous (piecewise constant) trajectories. In the modern literature, however, it is a standard convention to ignore the behavior of Lévy walks between the switching times, thus (mentally) substituting periods of motion with a constant velocity by the corresponding waiting time and an instantaneous transition (a jump). From this point of view one can look at Lévy walks as examples of general CTRWs with i.i.d. sequence of pairs $(X_i,W_i)$, where $X_i$ and $W_i$ are dependent. Namely, $X_i=V_i{W}_i$. The transition operator of the corresponding Markov chain gets the following form: with a finite measure $Q\left(ds\right)$ and a probability measure $P\left(dv\right)$.

In case of position depending Lévy walks, the transition operator takes a more general form with some stochastic kernels $Q(x;ds),P(x;dv)$. To simplify exposition we will not consider the extension when these kernels are time-dependent. As for general CTRWs with dependent distribution of jumps and waiting times one distinguishes the lagging Lévy walks $X_x\left(N_K\right)$ and the leading (or overshooting) walks $X_x(N_K+1)$.

As was mentioned, an appropriate scaling of CTRWs leads to a subordinated Markov process with averages evolving according to certain fractional (pseudo)differential equations. The general scheme for deriving such equations in case of i.i.d. sequences $(X_i,T_i)$ with dependent $X_i,T_i$ was given in [15] based essentially on the method of the Fourier and Laplace transform. The concrete version of this scheme for Lévy walks was performed in [16,17,18], which led to the so-called material fractional derivatives.

The general scheme for deriving fractional equations for scaling limits for CTRWs with position dependent jumps was developed in [19,20]. For position dependent jumps the methods of Fourier transform cannot be effectively applied, and completely different method had to be used. Modifications and extensions of this method were applied in [21,22] for equations with variable orders and kinetic equations. Here we apply this method to derive the equations for the scaling limit of Lévy walks for position dependent distribution of times and velocities (including variable orders of stability for the times of walks) leading to new general equations with fractional material derivatives.

In paper [23] it was suggested to look for various scaling limits for general CTRWs specified by the way the particles are considered to cross the boundary at the final jump (including leading and lagging processes as particular cases) and leading to equations with different multi-dimensional extensions of Caputo-Dzherbashian fractional derivatives. Here we argue that, specifically for Lévy walks, the most natural crossing rule is neither lagging nor leading process (where the last jump is included or not, respectively, in the final spatial position). In fact, as was already noted above, CTRWs are only approximations to Lévy walks, where particles are supposed not to jump, but to move with constant velocities between switching times. Thus the natural stopping time should be not before or after the final jump, but in the intermediate time, when the crossing of the boundary really occurs. This version of stopping was considered in some detail in [23] for general CTRWs and the corresponding fractional derivatives were derived under certain technical assumptions. Here we analyse in detail the equations governing the lagging, leading and intermediate limiting processes concretely for Lévy walks and construct the underlying limiting processes and Feller semigroups.

The distribution of waiting times $T_i$ in basic CTRWs (or moving times in Lévy walks) are assumed to have heavy tails with power decay at infinity and with unbounded expectation. In order to take into account possible position dependence, our main assumption will be that the time of moving T when started from a position x will have the power law with some positive function $\beta \left(x\right)$ bounded from above and from below. To make formulas more transparent we shall make (3) more precise (though this simplification is not very important). Namely, we assume that these distributions have continuous densities $Q_x\left(r\right)$ such that where $\beta \left(x\right)$ is some continuous function such that ${\beta }_1\le \beta \left(x\right)\le {\beta }_2$, with some constants $B>0$ and $0<{\beta }_1<{\beta }_2<1$ such that ${\beta }_1{B}^{{\beta }_1}>1$ (the latter condition ensures that ${\int }_B^{\infty }{r}^{-1-\beta \left(x\right)}dr<1$ for all x).

In the usual CTRW scaling (see e.g. [2,16,20]) one scales the transition times of the Markov chains (1) or (2) by some small parameter $\tau$ and the times of walks by $\tau$ in the power that equals the inverse value of the stability index of their distributions. Thus the natural scaling of jumps depends on the tails of jump distributions. Here we study Lévy walks without scaling velocities, so that the spatial scaling of jumps arises exclusively from the scaling of moving periods. The corresponding scaled version ${(X^{\tau },W^{\tau })}_{x,w}\left(t\right)$ of the Markov chain (2) can be defined by its governing transition operator

Here and the scaled inverse process is $N_T^{\tau }=sup\{t:W_{x,w}^{\tau }\left(t\right)\le T\}$.

The lagging and leading scaled CTRWs are then defined as the processes respectively.

As was mentioned, both these processes effectively neglect the fact that Lévy walks, strictly speaking, are not jump processes. Away from the boundary, this discrepancy is not essential, but for jumps crossing the boundary it becomes essential. In reality, if at a time $\tau (i-1)$ the process was at $(X_{x,w}^{\tau }\left(\tau (i-1)\right),W_{x,w}^{\tau }\left(\tau (i-1)\right))$ and a new moving period, say $W_i$, was revealed, which turned out to be final (that is, crossing the boundary$\{w=T\}$, or equivalently, such that $W_{x,w}^{\tau }\left(\tau i\right)>T$), then simultaneously a new velocity, say $V_i$, was revealed. Then the process starts moving from $X_{x,w}^{\tau }\left(\tau (i-1)\right)=X_{x,w}^{\tau ;T,lag}$ with the velocity $V_i$, until it reaches the spatial position of the intermediate exit point at the time $N_T^{int}=\tau (i-1)+\tau s$, when the process W reaches the boundary $\{w=T\}$, that is, where $s\in (0,1)$ solves the equation

In other words, the position $(X_{x,w}^{\tau }\left(N_T^{int}\right),W_{x,w}^{\tau }\left(N_T^{int}\right))$ in ${\mathbf{R}}^d\times \mathbf{R}$ is the point, where the straight line connecting $(X_{x,w}^{\tau }\left(\tau (i-1)\right),W_{x,w}^{\tau }\left(\tau (i-1)\right))$ and $(X_{x,w}^{\tau }\left(\tau i\right),W_{x,w}^{\tau }\left(\tau i\right))$ crosses the hyperplane $\{w=T\}$:

The transition operators for lagging, leading and intermediate stopped processes are the modifications of the transitions (5) for the initial Markov chain that take into account the chosen way of stopping at the boundary. Namely, they are defined for $w\le T$ and take the following forms:

We shall use the unified notation $U_{\tau }^{T,*}$ for these three operators, with * denoting either lag, or lead, or int. It is seen that all $U_{\tau }^{T,*}$ do not take the process away from any band ${\mathbf{R}}^d\times [a,T]$ for any $a<T$. Choosing $a=0$ without loss of much generality we shall consider, from now on, the corresponding processes as taking values in the band ${\mathbf{R}}^d\times [0,T]$ with some fixed T. Notice also that so that the corresponding processes are automatically stopped when reaching the boundary $\{w=T\}$.

The objective of the present paper is to analyse the limits of the discrete Markov chains ${\left(U_{\tau }^{T,*}\right)}^{[t/\tau ]}$, as $\tau \to 0$, and the corresponding subordinated processes (first coordinate subordinated by the inverse process to the second coordinate), to derive multidimensional Caputo-Dzherbashian-type and Riemann-Liouville-type fractional equations (with fractional material derivatives of variable order) that govern the evolution of the limiting processes, and to show well-posedness of these equations in natural functional spaces. Our approach is to avoid any difficulties arising from non-Markovian processes by first incorporating stopping rules in the Markov processes (5), identify the corresponding limiting generators and then look at the distributions of the final non-Markovian process as at the stationary distributions of the corresponding stopped Markov process.

We shall also discuss certain modifications of the model, namely when there can be additional waiting times between the walks, and when the walks are performed with parameter depending velocities.

The paper is organised as follows. In Section 2 we obtain preliminary results on the limiting Feller process for the scaled Markov chains (5) governing the sizes and times of jumps. Section 3 is the main one. It is devoted to the presentation of our main results concerning fractional equations that govern the scaling limits of the Markov chains (9) -(11), as well as the related process killed at an attempt to cross the boundary $\{w=T\}$. The latter turns out to be described by a multidimensional version (19) of the Riemann-Liouville fractional operator with material derivative, while the former are given by 3 different multidimensional versions (23) -(26) of the Caputo-Dzherbashian fractional derivative. We derive the corresponding equations and obtain their well-posedness. These results are contained in Theorems 3.4 and 3.5. Last SubSection 3.5 touches briefly on some modifications arising when using constant accelerations (rather than constant velocities) between the switching times, or even more general model with parameter depending velocities. The proofs of the main results are also essentially given in Section 3, up to some technical results, Theorems 3.1 - 3.3. The latter theorems, on the existence of the limiting Feller processes interrupted on an attempt to cross the boundary, are only formulated in SubSection 3.3. In order not to interrupt the main arguments, their proof is postponed and given in special Section 4. Section 5 is devoted to modifications arising from the inclusion of additional waiting times.

We conclude the introduction with certain notations that will be used in the paper without further reminder.

For a set $\Omega$, which is either ${\mathbf{R}}^d$ or the band ${\mathbf{R}}^d\times [0,T]$ in ${\mathbf{R}}^{d+1}$ (with some $T>0$) let $C\left(\Omega \right)$ denote the spaces of bounded continuous functions on $\Omega$, equipped with the standard sup-norm $\parallel .\parallel$. For $k\in \mathbf{N}$, let $C^k\left(\Omega \right)$ denote the spaces of k times continuously differentiable functions on $\Omega$ (in the case of a band, the corresponding one-sided derivatives are meant on the boundary) with bounded derivatives equipped with the standard norm where $\parallel f^{\left(m\right)}\parallel$ denotes the maximum of sup-norms of all partial derivatives of f of order m. Let $C_{\infty }\left(\Omega \right)$ denote the closed subspaces of $C\left(\Omega \right)$ consisting of functions vanishing at infinity, $C_{\infty }^k\left(\Omega \right)$ the closed subspace of $C^k\left(\Omega \right)\cap C_{\infty }\left(\Omega \right)$ consisting of functions such that all its derivatives up to order k belong to $C_{\infty }\left(\Omega \right)$.

Integrating by parts, (14) can be rewritten as

Recalling the expression for the standard (right) fractional derivative of an order $\beta \in (0,1)$: one can say that the integral represents (up to a positive multiplier that we shall neglect) the fractional material derivative, where the material derivative (in the direction v) is defined as

Thus we can write down the generator L in the form where the r.h.s. is the averaged (over v) fractional material derivative.

Continuing the analogy, let us note that if $f\in C^1\left(\mathbf{R}\right)$ vanishes at $w=T$, then the right Riemann-Liouiville derivative can be defined as the restriction of $D_{-}^{\beta }$ to the space of functions vanishing for $w\ge T$, that is, as the operator

Similarly, operators (16) and (17), reduced to the space of smooth functions vanishing for $w\ge T$, take the forms and

Thus $L_{T-}$ represents a multidimensional analog of the standard Riemann-Liouville fractional derivative (of variable order in our case), so that the inverse of this operator (if well defined) will represent a multidimensional analog of the standard Riemann-Liouville fractional integral.

When looking for a probabilistic interpretation of fractional derivatives in [23], it was noted that the Riemann-Liouville derivative is obtained from the free (without boundary) fractional operator by restricting it to the space of functions vanishing identically beyond the boundary, which in terms of the underlying stochastic process means its killing on the attempt to cross the boundary. In its turn, the Caputo-Dzherbashian derivative is obtained from the free fractional operator by restricting it to the space of functions that are constant beyond the boundary, which in terms of the underlying stochastic process means its stopping on the attempt to cross the boundary. This interpretation of fractional derivatives led to the natural extension of the fractional derivatives not only to a two-sided case, but to a variety of multidimensional settings. However, while killing on the boundary has always clear meaning, stopping for a multidimensional jump-type process really depends on the way one projects the result of the final jump (that crosses the boundary) to the boundary, which leads to several different version of the Caputo-Dzherbashian fractional operators. Three of them are analysed in this paper.

Using Lemma 2.1 allows us to conclude that the limits of as $\tau \to 0$, exist for $F\in C^1({\mathbf{R}}^d\times [0,T])$ and equal, respectively,

One can rewrite these expressions in the following equivalent forms:

Integrating by parts yields which is exactly the averaged fractional material derivative (18).

We shall denote by $L^{T,*}$ the general operators with * denoting either lag, or lead, or int. Our main technical results, given below, concern the existence of well defined Feller processes in ${\mathbf{R}}^d\times [0,T]$ generated by $L^{T,*}$.

To work with these operators let us denote by $C_{\infty ,0}({\mathbf{R}}^d\times [0,T])$ (or sometimes shorter by $C_{\infty ,0}$) the subspace of $C_{\infty }({\mathbf{R}}^d\times [0,T])$ of functions vanishing at the boundary $\{t=T\}$, and by $C_{\infty ,0}^1({\mathbf{R}}^d\times [0,T])$ the subspace of $C_{\infty }^1({\mathbf{R}}^d\times [0,T])\cap C_{\infty ,0}({\mathbf{R}}^d\times [0,T])$ with all partial derivatives belonging to $C_{\infty ,0}({\mathbf{R}}^d\times [0,T])$.

The elementary properties of $L^{T,*}$ are collected in the following proposition, its proof being obtained by a direct inspection that we omit.

We are also interested in the versions of these processes killed on the boundary $\{w=T\}$. The semigroups arising from killed processes act on the space $C_{\infty ,0}({\mathbf{R}}^d\times [0,T])$. It is seen that in this space all three operators $L^{T,*}$ coincide. Let us denote them by $L^{T,kill}$:

As expected, this is nothing else, but the operator (19), which represents (up to a sign) a multi-dimensional analog of the Riemann-Liouville fractional operator.

The proof of all these results follow the same line of arguments. We shall give details for Theorem 3.1 in Section 4 and briefly comment on modifications arising in other cases.

The following result is a straightforward but important corollary of the theorems given above.

Finally, when working with $L^{T,lead}$, we shall need to use functions from the domain that are not differentiable up to the boundary. Let us introduce the following functional space $H_{\infty }({\mathbf{R}}^d\times [0,T])$, which is the subspace of $C_{\infty }({\mathbf{R}}^d\times [0,T])$ consisting of functions $F(x,w)$ such that the spatial gradient ${\nabla }_{x}F$ exists and belongs to $C_{\infty }({\mathbf{R}}^d\times [0,T])$ and, with respect to the second variable, F is locally Hölder in the sense that the function is well defined and belongs to $C_{\infty }\left({\mathbf{R}}^d\right)$. It is seen that for any $F\in H({\mathbf{R}}^d\times [0,T])$ formulas (20), (21), and (22) yield well defined functions from $C_{\infty }({\mathbf{R}}^d\times [0,T])$. Consequently, using the fact that the generator of any Feller semigroup is a closed operator, and approximating the functions from $H_{\infty }({\mathbf{R}}^d\times [0,T])$ by the functions from the corresponding invariant cores of $L^{T,*}$ (given by the above Theorems) we obtain the following fact.

Let us start with the analogs of the Riemann-Liouville fractional operators.