{'fraction_non_alphanumeric': 0.04876921137060958, 'fraction_numerical': 0.013337629768692997, 'mean_word_length': 4.466227599938735, 'pattern_counts': {'":': 3, '<': 39, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
V-Langevin Equations, Continuous Time Random Walks and Fractional Diffusion

The following question is addressed: under what conditions can a strange diffusive process, defined by a semi-dynamical V-Langevin equation or its associated Hybrid kinetic equation (HKE), be described by an equivalent purely stochastic process, defined by a Continuous Time Random Walk (CTRW) or by a Fractional Differential Equation (FDE)? More specifically, does there exist a class of V-Langevin equations with long-range (algebraic) velocity temporal correlation, that leads to a time-fractional superdiffusive process? The answer is always affirmative in one dimension. It is always negative in two dimensions: any algebraically decaying temporal velocity correlation (with a Gaussian spatial correlation) produces a normal diffusive process. General conditions relating the diffusive nature of the process to the temporal exponent of the Lagrangian velocity correlation (in Corrsin approximation) are derived.


Introduction
Strange Transport has been the object of intense studies in recent years (A very recent qualitative review is [1]). (We use the terminology "Strange transport" [2] rather than "Anomalous transport" which is customary in Dynamical Systems theory, but has a different meaning in Plasma transport theory). The concept first appeared in the theory of stochastic processes, especially the theory of Continuous Time Random Walks (CTRW) [3], [4]. Consider a disordered system (e.g. a turbulent fluid or plasma). The position x(t) at time t of one of its particles is determined by its interactions with the other particles and/or with external sources. In a strongly disordered system these interactions are modelled by a random field. The statistical description of the system thus involves the definition of an ensemble of realizations. One of the important quantities describing the system is the mean square deviation (MSD) of the random variable x(t). In many cases, this is represented asymptotically (t → ∞) by a simple increasing function of time: where δx(t) = x(t) − x(t) . The value of the "diffusion exponent" µ characterizes the diffusion regime of the system 1 : The subdiffusive and the superdiffusive regimes are the ones called "STRANGE". It will be seen in the forthcoming sections that the concept of "strange transport" involves much more than a simple statement about the behaviour of the MSD.
All present theories of strange transport are of stochastic nature. From the very voluminous literature we only cite here a few of the more comprehensive, physically oriented review papers and books: [3] - [13], where additional references will be found.
As stated in [2], transport theories (e.g., for fluids or plasmas) can be constructed on three levels. α) A purely statistical mechanical theory would be based on the kinetic equation for a set of interacting charged particles, combined with Maxwell's equations. This would be the most fundamental description, but becomes impossibly complicated in practice.
β) The next level of description would be a compromise, based on a "semidynamical" model of particles moving according to the laws of mechanics (Newton's equation), but under the action of a fluctuating field, representing the action of the turbulent environment. This leads to stochastic ordinary differential equations of the Langevin type. More specifically, we consider V-Langevin equations [2], i.e. equations for the position of a "particle": where the right hand side represents a fluctuating velocity field, defined by its statistical properties. We always consider a divergenceless velocity field: Associated with (3) is a kinetic-type first-order partial differential equation for the fluctuating particle distribution function f (x, t) (whose average is the density profile), called a hybrid kinetic equation (HKE) [2]. By definition, the characteristics of this equation are precisely the V-Langevin equations (3): Due to the property (4) this equation is equivalent to the following: This equation ensures the conservation of the number of particles, and is rightly interpreted as a kinetic equation.
γ) The last level of description is the level of continuous time random walk (CTRW) theories and fractional diffusion equations (FDE), in which the deterministic dynamical laws are completely given up, and replaced by a purely random process.
All these concepts will be defined and used in the forthcoming sections. These three levels are, of course, interrelated: each one of them should be justifiable as an approximation of the more fundamental one.
The strange transport theories associated with CTRW and with FDE have been very successful in recent years in modelling many peculiar behaviors observed, among other fields, in plasma and fusion physics. There is, however, still a serious gap in justifying these purely stochastic models on the basis of a molecular theory, i.e., theories of the type α), or at least β).
In [2] we treated several models exhibiting strange transport, showing that a stochastic CTRW can be associated with a semi-dynamical V-Langevin equation under certain limiting conditions. All these models (diffusion in a fluctuating electrostatic field or in a fluctuating magnetic field) were either diffusive, or subdiffusive. One of the purposes of the present work is to try to determine whether there exist V-Langevin equations leading to temporal superdiffusive behavior and associating them with a CTRW or with a fractional diffusion equation.
A first work following the same philosophy was done by West et al. [14]. In that paper a specially simplified model of a one-dimensional, two-state CTRW was used. This means that a "particle" performs a CTRW in one dimension, moving with a velocity that is constant in absolute value: ±V . The particle moves during a given (random) time, after which it suddenly jumps and reverses its velocity. The velocity autocorrelation function is given as an algebraically decaying function of time. In a first step of their work, the jumps are supposed to be instantaneous events. This model leads already to possible superdiffusion of time-fractional type. The authors are, however, looking for a superdiffusive regime of Lévy type (space-fractional diffusion). They then modify their model. The jump PDF and the waiting time PDF of the CTRW are now linked by assuming that the particles move during their "jump" with final velocity ±V : thus longer jumps are "penalized" by longer durations ("Lévy walk " rather than "Lévy flight", [5]).
It appears to us that the limitation to two states makes the model unadapted to describing the motion of a physical particle; it is, however, necessary for the definition of the Lévy walk in their model. In the next section we show that the assumption about the two states can be dropped, and the velocity can be considered as a general random function of time: v(t) [Secs. 4 -5].
The model can be further extended by considering a velocity field : v(x, t) (Secs. 6 and following). But this case leads to additional difficulties, because it requires the determination of the Lagrangian velocity correlation, a problem that is well known as a basic difficulty of turbulence theory. Here we shall consider only the crudest approximation for the latter quantity, but the application of more precise methods that were recently developed in this field can be envisaged.

CTRW and Fractional Calculus
In the forthcoming text, the abbreviation "PDF " denotes a "probability distribution function"; the symbol may equivalently be translated as "prob-ability density function".
In the Continuous Time Random Walk model one considers a particle which at time zero makes an instantaneous jump from position x 0 to x 0 + x, then waits at the new position during time t, then jumps again to a new position, where it waits, etc. Both x and t are random quantities. In all forthcoming calculations, we assume that the random processes are symmetrical (centred), i.e., x(t) = 0, ∀t; hence [see Eq. (1)], δx(t) = x(t). In its standard form, a CTRW is defined by two functions: f (x): the PDF of a jump defined by the vector x; Fourier transform: f (k), ψ(t): the waiting time PDF; Laplace transform: ψ(s).
We only consider in the present work "classical" CTRW's, for which f (x) and ψ(t) are independent quantities. This means that f (x) is independent of t, and ψ(t) is independent of x. Processes such as Lévy walks [5] or the vMSC model [15] are not discussed here.
A derived function that plays an important role is the memory kernel φ(s) (in Laplace representation): This standard CTRW has been solved exactly by [3] in 1965. Let n(x, t) be the PDF of finding the walker at position x at time t, given that it started certainly at the origin x = 0 at time t = 0 [i.e., n(x, 0) = δ(x) and n(x, t) → 0 as |x| → ∞]. Clearly, from this definition, n(x, t) is the propagator (or Green's function) of the CTRW process. It is given, in Fourier-Laplace representation, by This is the celebrated Montroll-Weiss equation. An equivalent form is: The latter equation is easily transformed into: This is the Fourier-Laplace transform of the integro-differential equation obeyed by the propagator: This Montroll-Shlesinger (MS) Master equation is non-local, both in space and in time. It shows that the rate of change at point x and time t is determined both by the past history, whose influence is measured by φ(t), and by the spatial environment characterized by f (x).
Although Eqs. (8) and (9) are explicit solutions for the propagator (or the resolvent, as the Fourier-Laplace transform of the propagator is usually called), the expression of the inverse Fourier-Laplace transform in terms of known functions of x and t is, in general, difficult or impossible.
It has been shown in a series of important recent works that this problem can be solved under rather general conditions in the so-called fluid limit, which is actually an asymptotic limit of large distances and times [ [5], [6], [13], [16], [2]]. (It should be noted that, in many interesting cases, the notion of "long" is ambiguous when there is no intrinsic time scale!) The fluid limit is defined as follows (using dimensionless variables k, s, x, t) [d is the dimensionality of the system]: In the same limit, the memory function is: This is the same as Eq. (17.12) of [2], describing subdiffusion (β < 1) or normal diffusion (β = 1) [see Sec . 3]. Later in the present work we will need an extension to the range 1 < β < 2. This extension is not trivial and requires a closer discussion, which is developed in Sec. 4. Thus, in the fluid limit, the CTRW is entirely defined by the two exponents α and β.
Its solution, which is the resolvent of the process, is: Eq. (15) is, in F-L form, a fractional differential equation (FDE), corresponding to: (17) A primer on fractional calculus, including the definition of these objects is given in the Appendix A. Thus, Eq. (15) can be rewritten as: which is the F-L transform of the FDE: with the initial condition: Thus, n(x, t) is the propagator (or Green's function, or fundamental solution) of the FDE (19).
Although the Fourier-Laplace inversion of Eq. (16) is difficult, an important property, i.e., the scaling of the solution in x − t space can be derived without performing explicitly the latter operation. The following argument is due to Mainardi [9], and is much more compact than the one presented in [17], [18] and [2]. It is based on the well-known properties of the Fourierand Laplace transformations: where the double arrow connects a pair of Fourier (Laplace) transforms. Applying this property to Eq. (16), we find, after a short calculation: Applying the relations (21) to the last form we find: Thus, finally: Let us assume now that the density profile has the scaling form: Under what condition is this form compatible with Eq. (23)? The left hand side of the latter yields: and the right hand side yields: These two expressions are equal, provided that γ = β/α. We thus conclude: This is the important scaling relation obeyed by the density profile. It tells us, in particular, that it depends on space only through the combination p = (x/t β/α ), called the similarity variable. From this result we immediately deduce the scaling of the MSD: With the substitution p = x/t β/α we find: identical to Eq. (1), with the following definition of the constant: M α,β = dp p 2 K α,β (p).
Eq. (26) provides us with a criterion of "strangeness" expressed in terms of the two parameters α, β: The dimensional argument leading to Eq. (26) does not tell us anything about the constant M α,β : we do not even know wheter it is finite or infinite! The determination of M α,β requires the explicit form of the solution through the function K α,β (p).
The first terms of the expansion of the Fourier-density profile are thus: The MSD can be determined directly from the Fourier-density profile by the well-known formula: When applied to the Mittag-Leffler function, this yields (in d dimensions): Two cases are to be distinguished: In this case, the first term has a finite positive value in k = 0; all the other terms in the series vanish in this point. Hence: This completes Eq. (26) The first term in the right hand side of Eq. (34) contains the factor |k| α−2 which diverges in k = 0. The remaining terms are either 0 or ∞ in this range of α. Thus, the constant in Eq. (26) is: This means that the density profiles of the CTRW in the present range of α have long tails, whatever the value of β.Restricting the concept of superdiffusion to the sole discussion of Eq. (26) would put all the processes with 0 < α < 2 on the same rank. The distinction among them requires a deeper discussion of the density profile. These questions will be further discussed in the next Section.

Classification of Strange Transport Regimes
The very clear papers of Mainardi and coworkers [8], [9], [10] contain an exhaustive study of the Green's function of the FDE (19) in one dimension, d = 1. In the forthcoming work, we take full advantage of the fact that the general solution n α,β (x, t) (for all α, β) of Eq. (19) has been obtained by these authors.. This solution was also obtained by [5] in a different form (in terms of Fox functions). We shall first discuss some particular cases (for d = 1).
This is the classical case of normal diffusion. Eq. (15) becomes in this case: which is simply the F-L transform of the ordinary diffusion equation (with a diffusion coefficient D = 1): It is well known that in this case the MSD is: in agreement with Eq. (26). The solution of Eq. (36) is: The propagator associated with Eq. (37) is the Gaussian wave packet: Note that Eq. (40) is a special case of the scaling relation (25): The Gaussian scaling function F (p) in the present case is a consequence of the Central Limit theorem of probability theory. "Most" physical systems behave in this normal diffusive way. But it is the strange behavior, which deviates from this Gaussian form, that will be the main object of our study.
This case was studied in detail in [2], and more completely in [8] and [9]. It was called the Standard Long Tail CTRW (SLT-CTRW) in [2], and time-fractional subdiffusion in [9]. The main result in this connection may be stated as follows. Eq. (26) shows that the SLT-CTRW is a subdiffusive process (because 2β < α), for all values of β in the range (0, 1). The propagator is non-Gaussian, and has a long tail of stretched exponential form: The values of the constants A α,β and b β [9] are irrelevant here. Eq. (42) has again a scaling form, with a (rather complicated) finite constant A 2,β . This function has a single maximum in x = 0, and has a finite MSD.
This is the case of space-fractional diffusion of Mainardi [9] [see also Sec. 12.3 of [2]]. The corresponding resolvent is: Its inverse Laplace transform is: This is the Fourier transform of a Lévy stable distribution ([2], Appendix D). It is well known that the inverse Fourier transform of these distributions is in general not expressible in terms of simple functions, except for a few values of α. The general propagator (for arbitrary α) can, however, be expressed as a Mainardi function, defined by a convergent power series [see Eq. (77) below]. The most striking property of the Lévy distributions is their "fat" long tail, i.e., an algebraic power law decay for p ≫ 1: This law of decay is much slower than the stretched exponential decay of Eq. (42). As a result, the MSD of all the Lévy distributions is infinite, as shown in Sec. 2. The three cases discussed above are illustrated in the following figure by the corresponding Mainardi functions. In these figures t = 1, hence p = x.
We thus see that there are two types of long tails: • The stretched exponential tail, characteristic, typically, of the subdiffusive SLT • The algebraic tail, characteristic of the superdiffusive Lévy process.
Note that both are "fat" tails, i.e. the decay at large x is slower than the corresponding diffusive Gaussian. The second type is the most slowly decaying one. Warning: it should not be concluded that all subdiffusive processes have exponential tails, whereas all superdiffusive ones have algebraic tails!! (see below).
This case corresponds to time-fractional superdiffusion. This very interesting case will be treated in detail in Sec. 5.
This is the generic, mixed spatial and temporal fractional diffusion. In this case there is a competition between super-and subdiffusion, expressed by the ratio 2β/α [see Eq. (26)]. In the work of del-Castillo Negrete et al [16], a model with α = 3/4, β = 0.5 has been successfully compared with numerical simulations of a set of plasmadynamical equations. In a recent paper [21], additional values of α < 2 were explored.
In the present work we do not consider the case of space-fractional (Lévy) superdiffusion. We thus fix in all subsequent sections: 4 Random Time-dependent Velocity We consider a particle moving in one dimension; its position X(T ) at time T is determined "semi-dynamically" by a V-Langevin equation ( [2], Ch. 11): 2 where the random velocity V (T ) (depending solely on time) is defined statistically by specifying the first two moments of the ensemble of its realizations, supposed to be Gaussian, homogeneous and stationary: The function T (T ) need not yet be specified explicitly, except by requiring T (0) = 1. V 0 , the initial velocity rms is supposed to be a given parameter.
Let f (X, T ) be the (fluctuating) distribution function of the particle. It is represented as usual as: where the ensemble-average n(X, T ) ≡ f (X, T ) is called the density profile, and δf (X, T ) is the fluctuation. The distribution function obeys the Hybrid Kinetic Equation (HKE): The characteristic equation of this first-order partial differential equation coincides with the V-Langevin equation (47). This equation is split, as usual [2], into two components by the averaging operation: The latter equation is solved by using the propagator: Substituting its solution into (51), we find: This is the exact solution of the problem. It is expressed by an integrodifferential equation that is non-local both in time (non-Markovian) and in space (non-local). We do not repeat here the arguments justifying the local approximation, by which the fluctuating displacements in the arguments of n and of δf are neglected, and by which the second term of the right hand side, containing the initial fluctuation, is also neglected [2]. The result is an (approximate) local, but non-markovian diffusion equation: Using Eq. (48), this becomes: Consider, on the other hand, a particle performing a one-dimensional CTRW characterized, in the fluid limit, by a Gaussian jump PDF: where σ is a characteristic length and κ is the dimensional wave vector. The form of the waiting time PDF is not specified at this stage.
The density profile 3 (in F-L representation) obeys Eq. (10) which, written with a dimensional Laplace variable S, combined with (57) reduces to: Its inverse F-L transform is: This equation constitutes the main result of the first part of our work. Comparing it with Eq. (56), we see that the density profile produced (in the local approximation) by the semi-dynamical V-Langevin equation with a time-dependent V (t), and the density profile produced by a CTRW with α = 2, -or, more specifically, by the fractional diffusion equation (58)obey formally identical equations of evolution. This statement will be made more precise in the discussion of Sec. 5. The quantitative correspondence is obtained by identifying the velocity correlation function of the V-Langevin process with the memory function: i.e., in F-L representation We note that Eqs. (56) and (58) were also obtained in Ref. [14] in a quite different manner, which is only valid for the 2-state CTRW. The relation (61) is much more general and was obtained in a simpler direct way.

Algebraic Velocity Correlation
We now make the problem explicit, by adopting the same form of the velocity correlation as in Ref. [14]: where A is a quantity having dimension [length] 2 [time] −2+γ . We then associate with the V-Langevin equation (for d = 1) (47) [or the Hybrid kinetic equation (50)] a CTRW with memory function given, in the fluid limit, by Eq. (60), leading to the fractional diffusion equation (FDE) (58): This relation agrees with Eq. (33) of [14]. It is trivially equivalent to: The F-L density profile is obtained from Eq. (63) in the form: Comparing this result with Eq. (16), we immediately see that the two expressions become identical by taking the following value for the waiting time exponent β: Thus, the process corresponds to time-fractional superdiffusion as defined in Sec. 3: At this point there appears a subtle difficulty (the following argument results from a remark of J. Klafter [22]). The FDE (63) is formally derived as the fluid limit of a CTRW with ψ(S) = 1 − S β , and β in the range (66): this is an extension of the range defined in Eq. (13).
The knowledge of the Laplace transform ψ(S) allows us to calculate the average waiting time between two jumps in the CTRW, by using a well-known formula: For 0 < β < 1 the average waiting time of the CTRW is infinite. This denotes a waiting time PDF with a long tail, leading to strange transport, i.e., subdiffusion (Sec. 3 and [2]). For β = 1 the average waiting time is finite, corresponding to normal diffusion.
For 1 < β < 2, the average waiting time is zero, thus leading to a paradoxical situation. Recalling that T is a semi-definite positive variable, it follows that ψ(T ) ∼ δ(T ). Thus, any random walker would remain trapped forever at its starting point, in clear contradiction with superdiffusion, as predicted by the diffusion exponent µ = 2β/α = β > 1 (Sec. 3). The only other possibility of obtaining a zero average waiting time would be to have ψ(T ) negative in part of the range (0, ∞), but then ψ(T ) does not qualify for a probability distribution function.
The conclusion of this argument is the following. The function n(κ, S) defined by Eq. (65) is, indeed, the resolvent of the fractional differential equation (FDE) (63), but the latter cannot be derived as the fluid limit of a CTRW.
We are now faced with a new subtlety 4 . The function n(κ, S) is the solution of the linear algebraic equation (64). But, is the latter equation the F-L transform of a fractional diffusion equation? One is tempted to associate Eq. (64) with the FDE in X − T representation: with the initial condition: This fractional equation is an "intermediate" between a diffusion equation (γ = 1) and a wave equation (γ = 0) [23], [8], [10]. But, using the general formula for the Laplace transform of the Caputo fractional derivative [7], [10] and the Appendix of the present paper, the F-L transform of Eq. (69) is: (71) It thus appears that the general solution of the Cauchy problem for Eq. (69) requires two initial conditions (like the wave equation): n(X, 0) and n T (X, 0) ≡ [∂n(X, T )/∂T ] T =0 . As a result, the general solution is expressed in terms of two propagators; indeed, the solution of Eq. (71) is: . We refer to Ref. [10] for a complete solution and discussion of this problem. Here we shall only consider the simpler special case (also considered in [9]): In this case, Eq. (71) reduces to Eq. (64), and Eq. (65) represents the F-L propagator of Eq. (69) equipped with the initial conditions (73). The explicit solution for the propagator of (69) (i.e., the F-L inversion of Eq. (65) was obtained analytically by Mainardi [23], [9] (and also by Metzler and Klafter [5] in terms of Fox H-functions): This has precisely the scaling form (25), with the similarity variable: The special function M 1−(γ/2) (p) is called the Mainardi function: it is defined for any order ν ∈ (0, 1) as: This is a special case of a Wright function, defined as follows ( [20]): The series (77) is convergent, but the convergence is very slow. In practice, for the numerical calculations displayed below, we take 100 terms in the series. Mainardi et al 2001 [9] also showed that the asymptotic behavior of the propagator is given by a stretched exponential, which for the value β = 2 − γ is: with: The important point here is that the present density profile decays asymptotically according to a stretched exponential law, with exponent 2/γ > 2. This implies that the propagator decays asymptotically faster than the Gaussian (contrary to the subdiffusive case 0 < β < 1, discussed in Sec.3. Mainardi calls this feature a "thin tail ". It thus appears clearly that superdiffusion is not necessarily associated with a long ("fat") tail of the propagator.
Another important property of the propagator (75) [i.e., of the Mainardi function M ν (p) for 1/2 < ν < 1, i.e., 1 < β < 2] is its possessing for T > 0 two symmetrical maxima that move away from the origin, while becoming wider. As pointed out above, the fractional equation (69), written in a simpler notation (∂ β /∂T β ≡ C 0 D β T ) as: is a kind of interpolation between a diffusion equation (β = 1) and a wave equation (β = 2). Thus, the solution represents a wave propagating in both directions away from the origin, and damping out on the way 5 . We illustrate these properties by the plot of Fig. 3, corresponding to γ = 0.5, hence β = 1.5. For T → 0, the propagator tends toward n(X, 0) = δ(X), as it should (70). We show in Fig. 4 the tail of the propagator (for T = 1) compared to the tail of the corresponding Gaussian. This situation is to be contrasted with Figs. 1 and 2: both the time-fractional subdiffusive propagator and the space-fractional superdiffusive (Lévy) propagator have "fat" tails.
Figs. 3 and 4 allow us to understand how superdiffusion is possible in spite of the thin exponential tail. The dispersion of the matter is produced not only by the broadening of the initial distribution (as in ordinary diffusion) but also by the symmetric wave-like outward motion of the maximum.
The results of the present section are not really new: All these results were also obtained in Ref. [14], using, however, the 2-state model. We showed here that the results are easily extended to an arbitrary fluctuating velocity depending solely on time. Our main result can be summarized as follows.
The motion of a particle in one dimension, with a random velocity depending solely on time, obeying a semi-dynamical V-Langevin equation (V-LE) (47), can be equivalently described by a time-fractional differential equation (FDE). There is, however, no CTRW whose fluid limit is the FDE (74). The V-LE is characterized by a velocity correlation depending (for all positive times!) algebraically on time, T (T ) ∼ T −γ , with an exponent γ ∈ (0, 1). The equivalent FDE is characterized by a spatial exponent α = 2 and a waiting time exponent β = 2 − γ; thus 1 < β < 2. This FDE leads to time-fractional superdiffusion, with diffusion exponent µ = 2 − γ > 1. The density profile has two maxima that are moving symmetrically away from each other, while broadening with increasing time. The asymptotic decay is exponential, decreasing faster than a Gaussian in space ("thin tail").
The model V-LE considered here, though allowing an exact analytical solution, has a physical drawback. The algebraic law of the velocity correlation T (T ) = T −γ cannot be valid for short times. For an obvious physical reason, lim T →0 T (T ) = 1 [as required by Eq. (48)], instead of diverging as in the present model. The model should certainly be improved by using a correlation function that satisfies this constraint.
We now generalize our model by considering, instead of a velocity depending on a single time variable, a velocity field in two dimensions, for which the velocity depends on both time and space V (X, T ).

Random Velocity Field with Algebraic Timecorrelation
We now return to Eq. (3), and consider a true velocity field, i.e., a function both of space and time. We immediately note that this generalization introduces nothing new in one dimension. We recall the important condition of zero divergence (3) which ensures the correct physical status of the hybrid kinetic equation, i.e., the equivalence of Eqs. (4) and (5). Clearly, the constraint (3) reduces, in one dimension, to ∂ X V (X, T ) = 0, hence the velocity can only depend on time, and we are brought back to the situation of Sec. 4. In order to study a non-trivial problem, we must consider d ≥ 2. Specifically, we shall consider in the rest of the present work two-dimensional systems, d = 2. As will be shown below, such systems are actually of much greater physical interest than one-dimensional ones. We thus consider a particle whose position is specified by the vector X = (X, Y ) and the velocity V = (V x , V y ). The corresponding (vectorial) V-Langevin equation is: As a result of the previous discussion, we may associate with this V-Langevin equation a hybrid kinetic equation, provided that the velocity field has zero divergence: Contrary to the one-dimensional case, this condition can be realized with a non-trivial velocity field. The associated hybrid kinetic equation possesses both required properties: it conserves the normalization of f and its characteristic equations are Eqs. (82). Thus, the two following forms are equivalent: We note that the condition of zero divergence of a two-dimensional velocity field is realized, in particular, by the intensely studied model of a charged particle moving perpendicularly to a strong constant magnetic field B, in presence of a fluctuating electrostatic field (see, e.g. [2], where many other references are given]. The coordinates of its guiding centre are animated by the well-known electrostatic drift velocity (in Gaussian units): where Φ(X, T ) is the (fluctuating) scalar potential. The statistical properties of the velocity field are thus derived from the statistical definition of the potential. We assume that the latter is a Gaussian centred, homogeneous, isotropic and stationary random process, defined by a vanishing first moment and by the following factorized form of the second moment: The second relation defines the Eulerian potential correlation, i.e., the correlation between the values of the potential at the origin and at a fixed point X at time T . In the present work we consider only the case where the spatial part of the correlation has a Gaussian form: where λ is the so-called correlation length. The temporal part T (T ) will not yet be specified explicitly. It should, however, satisfy the condition: T (0) = 1, and is supposed to depend only on the absolute value of T .
It is convenient to introduce from here on dimensionless quantities (denoted by lower case letters), defined as follows: Here ε is a measure of the intensity of the potential fluctuations, e.g., the root mean square value of the potential fluctuations. T c is an intrinsic characteristic time, to be defined later [Sec. 8]. The V-Langevin equations of motion (86) reduce to: The dimensionless constant K is called the Kubo number: Given the form of Eqs. (90) it is even more convenient to use a further reduction of the time variable: The V-Langevin equations then reduce to: The relevant dimensionless Eulerian correlations are now defined as follows: with (88): A first, well-known general result is obtained from the V-Langevin equation (93): the mean square deviation (MSD) in the x-direction is expressed as follows: Here L xx (θ) is the Lagrangian velocity correlation, calculated along the trajectory x(θ) of the particle, i.e., using the solution of Eq. (93): The MSD is related to the dimensionless running diffusion coefficient D x (θ) in the x-direction by the well-known Einstein relation: 6 We now consider the evolution of the density profile n(x, θ) = f (x,θ . The distribution function is governed by the hybrid kinetic equation (84), or equivalently, (85). The latter is treated by a straightforward generalization of the reasoning leading from Eq. (50) to (55); it will not be repeated here. The final result differs from the latter equation in that the simple one- Clearly, the relation (96) can be obtained directly from Eq. (99).

The Corrsin Approximation
The study of transport is now reduced to the (difficult) determination of the Lagrangian velocity correlation. More precisely, we want to find a relation between the (given) Eulerian velocity correlation and the Lagrangian one.
In the present work we will use only the simplest approximation method for this problem: the Corrsin approximation [24], [2]. We rewrite the exact Lagrangian velocity correlation as follows: In the Corrsin approximation, the average in the integrand is assumed to be factorized as follows: The first factor in the integrand is recognized as the Eulerian correlation tensor, assumed to be of the form (87), thus: The last factor is evaluated as follows, in the second cumulant approximation: Thus: We now use the explicit forms (95) for the Eulerian correlation. It is immediately seen by symmetry that the non-diagonal components of the Lagrangian tensor vanish, and the two diagonal ones are equal to each other. with: The integration over x and y is elementary. The result is combined with Eq. (96) to yield: We thus ended with an integral equation for the Lagrangian velocity correlation in the Corrsin approximation [26], [2]. This equation will be solved in Secs. 9 and 10.
The diagonal character of the Lagrangian velocity correlation introduces a significant simplification in the equation for the density profile, which reduces to the following non-Markovian diffusion equation: The diffusion coefficients in the x-and y-directions are equal to each other. We may therefore drop the subscripts in Eq. (98) and write simply: Whenever this integral converges, its limit for θ → ∞ defines the ordinary diffusion constant: D = D(∞). Whenever D is finite and positive, the regime is normal diffusive.

Algebraic Time Correlation: General Properties
We now make the problem more explicit. We consider a class of V-Langevin equations (and associated hybrid kinetic equations) for which the temporal part of the Eulerian velocity correlation, T (θ/K), is a long-tailed function that decays algebraically for long times. It is of the general form (62), but we take now some additional precautions to ensure a correct behavior at short times. We must, in particular, avoid the divergence at T → 0. Physically, any correct (dimensional) correlation function should have the property: T (T = 0) = 1. We thus define a characteristic cut-off time T c as the value of T at which the right hand side of Eq. (62) equals 1. For all times 0 ≤ T ≤ T c we put T (T ) = 1. Note that a correlation of the form (62) has no intrinsic time scale; but the adoption of a cut-off introduces a characteristic time T c . The latter allows us to complete the definitions of the dimensionless time t = T /T c (89), of the Kubo number (91) and of the dimensionless variable θ (92). We thus adopt the final form of the dimensionless temporal Eulerian velocity correlation as: For convenience, γ will be called the Eulerian exponent; by definition it will be limited to the range (0, 1). Anticipating the results of the next Section, it will be shown that the Lagrangian correlation function L(θ) has the same general algebraic form, with different values of the exponent and of the coefficients: In order to ensure a correct connection of the two parts, we must take: Λ will be called the Lagrangian exponent: its value depends on the Eulerian exponent γ and on K. It will be determined in the next sections.
Before going into specific calculations, we wish to derive an important property relating the nature of the macroscopic transport to the value of the Lagrangian exponent. Combining the definition of the dimensionless running diffusion coefficient (108) with the form (110) we obtain: The integrals are easily evaluated: The more interesting long-time part is rewritten as the sum of a constant and of a time-dependent term; using also Eq. (111) we find: We now see the existence of two sharply separated ranges of the Lagrangian correlation exponent Λ.
A) 0 < Λ < 1. In this case the second term in the right hand side of (114) is a monotonically increasing function of time. Hence, for sufficiently long times this term strongly dominates the first constant term, and the running diffusion coefficient behaves asymptotically as: The regime is thus (time-fractional) superdiffusive.
The behavior is radically different. The second term in the right hand side of Eq. (114) is now a monotonically decreasing function of time. For long times this term decreases to zero. Hence the running diffusion coefficient tends toward a positive constant: The regime is now normal diffusive for all values of Λ > 1. The asymptotic D is the ordinary diffusion constant.
Clearly, the value Λ = 1 is a very special one: as Λ reaches 1 from below, the regime changes abruptly from a superdiffusive to a normal diffusive one, and remains so for all values of Λ > 1. This phenomenon is illustrated in Fig. 5. The difference between the normal diffusive regime obtained from the kinetic equation and the subdiffusive regime predicted by the CTRW for the same value of Λ > 1 appears clearly in Fig. 6.
In order to obtain a clearer picture of what happens at this point, we study the density profile. More specifically, we investigate whether its evolution can be described by an equivalent FDE, as was found in Sec. 2.
As shown in Sec. 4, there exists a formal equivalence between, on one hand, the non-Markovian diffusion equation (107) or, in Laplace representation: The Laplace transform of the Lagrangian correlation (110) is, for long times: Comparing with Eq. (14) we find that the temporal exponent of the "formally equivalent" FDE is: Substituting this value into Eq. (16) we find a result similar to (65), with Λ replacing γ. The resolvent given by the Montroll-Weiss equation is: The difference with (65) is that we should not limit consideration to Λ < 1. We consider, instead the various cases discussed above.
The situation is the same as in Sec. 5. There is indeed a complete identity between the kinetic result and the FDE with 2 − Λ, i.e., 1 < β < 2. As shown, however, in Sec. 5, there is no equivalent CTRW in this case. The regime is one of time-fractional superdiffusion. The density profile is provided explicitly by the Mainardi function (75) 7 . It has two moving maxima, as represented in Fig. 3 and decays asymptotically as a stretched exponential. The diffusion exponent µ = 2 − Λ. All this is in agreement with the kinetic result. B) Λ = 1. The resolvent (122) of the FDE reduces to: This is the Fourier-Laplace transform of a Gaussian packet, characteristic of a normal diffusive regime. On the other hand, the result (116) shows that the running diffusion coefficient obtained from the kinetic equation grows, slowly but indefinitely in time. The FDE is no longer equivalent to kinetic theory: Λ = 1 appears as a bifurcation point.
The propagator (122) corresponds now to a FDE with α = 2, 0 < β < 1. It describes a time-fractional subdiffusion, which can be derived as the fluid limit of a CTRW. The corresponding propagator has a single maximum at the origin, and decays according to a stretched exponential [see Eq. (42)]. This regime has been studied in detail in [2]. On the other hand, the kinetic theory predicts a normal diffusive regime throughout this range of Λ [Eqs. (114), (117)]. There appears now a divorce between CTRW, FDE and kinetic theory. Λ = 1 has indeed the property of a bifurcation point: as the parameter Λ increases from 0; the propagator first follows the "FDE branch". At Λ = 1 it leaves this branch and follows the completely different "diffusive branch".
The kinetic equation for the propagator in this regime is obtained from Eqs. (99) and (110): For long enough time, the equation can be Markovianized by neglecting the retardation in the right hand side. This approximation is justified by our knowledge that a constant diffusion coefficient exists. The final asymptotic equation is an ordinary diffusion equation: with D given by Eq. (117). Nevertheless, for any finite time this Markovianization is not reliable, because of the non-exponential Lagrangian correlation, hence the density profile is in general non-Gaussian. D) Λ > 2. In this range the divorce between kinetic theory and FDE is complete. The latter looses its meaning. Indeed, it would correspond to β < 0, which would describe a waiting time distribution that is an increasing function of time. Alternatively, its Laplace transform diverges: ψ(s) ∼ s −|β| → ∞ as s → 0. On the other hand, the Markovianization of the kinetic equation (124) is more and more justified as L(θ) decreases faster for large θ.

Lagrangian Velocity Correlation: First Approximation
We now determine the Lagrangian velocity correlation as a function of time and of the parameters γ and K of the input Eulerian velocity correlation, by solving the integral equation (106). The procedure used is one of successive numerical iterations L (n) (θ), starting with the natural Eulerian trial function L (0) (θ) = T (θ/K). We checked the convergence of the procedure. As can be seen in Fig. 7, the second iteration provides already an excellent approximation.
The first very important qualitative result is the following. An input Eulerian correlation (109) with an algebraic tail produces a Lagrangian correlation having the same type of algebraic tail. This is clearly seen in a log-log representation of the result, as shown in Fig. 8. We thus confirm the anticipated form (111). Our main task will be the determination of the Lagrangian exponent Λ, of the corresponding "intensity" h and of the cut-off time θ L as functions of the input Eulerian exponent γ and of the Kubo number K.
A first glance at Fig. 7 shows that, even for relatively small values of K, the Lagrangian correlation deviates very quickly from the Eulerian one: this is clearly seen in Fig. 8. Thus, the nonlinearity, measured by K, is far from producing just a small perturbation.
The first iteration is not quantitatively accurate (Fig. 7), but it has the advantage that it can be evaluated analytically. It is therefore interesting   to study explicitly some of its properties which provide a global qualitative picture. The explicit form is (106): The integral in the denominator, combined with the definition (109), can be written in terms of two other integrals as follows: where: Combining these partial results, we finally obtain: The behavior of this function for short times is very simple: starting from 1, it decreases parabolically: As θ reaches K, we find, of course, that the two branches (130) and (131) join each other at the same point. But the derivative at θ = K is discontinuous. Indeed, consider the two integrals (128) and (129); one easily finds: Thus, L (1) (θ) exhibits a kink at θ = K, as can be seen in Fig. 9. The asymptotic behavior for large θ is more interesting. Here we must distinguish two ranges of the Eulerian exponent γ. a) γ < 1. In this case the term proportional to θ 2−γ dominates all others in the denominator of Eq. (131) for large θ. We thus find: We thus find, indeed, for very small γ, an algebraic tail of the form (110): with: Next, we note that the Lagrangian exponent in this approximation is: Although this value is not accurate (as can be seen from Fig. 8), it has the important property (which will be confirmed) that it only depends on γ, not on K (for very small γ). b) γ > 1. In this case the dominant term at large θ in the denominator of Eq. (131) is the term proportional to θ. The asymptotic behavior is thus different: where k (1) (K, γ) is found from this equation. The important point is that the Lagrangian exponent takes a different dependence on γ for larger values of this exponent: Thus, Λ (1) (γ), which was a decreasing function, starts increasing with γ beyond γ = 1. We may calculate analytically the exact dependence of Λ (1) (γ) by evaluating the slope of the log-log graph of L (1) (θ) (see Fig. 8): The result is shown in Fig. 10. Comparing this type of figures for different values of K, one finds a very weak dependence of Λ on K for γ close to 1.
The important consequence of this new "bifurcation" is the fact that, for all values of γ we find: The first iteration thus predicts, in accordance with the discussion at the end of Sec. 8, that the first iterate of the Lagrangian correlation leads to a diffusive regime for all values of γ. More than that: the Lagrangian exponent is in the range D), where there is no equivalence between kinetic theory and CTRW or FDE. This statement, however, has to be confirmed by the final solution. 10 Lagrangian Velocity Correlation. Final Numerical Calculation.
As shown in the previous section, the second iterate of Eq. (106) yields already an excellent approximation, while not requiring long numerical calculations (we therefore omit the superscript (2) ).The short-time (θ < K) part can even be calculated analytically. We thus have: Eq. (131) is substituted in the integral, and the expression is evaluated numerically.
For short times, the Lagrangian correlation behaves as: i.e., exactly as in the first approximation (132). As can be seen in Fig.  11, the kink which appeared at θ = K in the first approximation is still present in the final expression. The long-time behavior is, of course, more interesting. The first important point, that is seen in Fig. 8, is that the Lagrangian velocity correlation behaves asymptotically as a power law, (110), thus confirming our previous anticipation. We now study its parameters.
The Lagrangian exponent Λ is in all cases smaller than the one given by the first iterate, but much larger than the Eulerian exponent γ, as seen in Fig. 8. In Fig. 12 we plotted log L(θ) vs log θ for given γ = 0.1 and three values of K = 0.1, 0.5, 0.7. For all these relatively small values of K 8 the exponent Λ is strictly the same: Λ(γ = 0.1) = 2.10, ∀K.
Next, we measured the dependence of Λ on the Eulerian exponent. The result is a very simple linear dependence, shown in Fig. 13. This dependence is very accurately fitted by the formula: We note already here the following very important point: the Lagrangian exponent Λ(γ) > 2, ∀γ.   We now consider the numerator h of Eq. (110), determined from the logarithmic representation of L(θ) (Fig. 12). We first note (Fig. 14) that, for given K, it is a very slowly varying function of γ: we therefore consider it approximately independent of γ. The dependence on K appears to be quadratic, for relatively small values of K. It is quite well fitted by the following empiric formula, valid for γ 0.5 (Fig. 15): Finally, the value of the cut-off time is obtained by setting L(θ L ) = 1.
To sum up, we found a reasonable analytical representation for the Lagrangian velocity correlation in the Corrsin approximation. A simple feature of this representation is the fact that (see Fig. 13 and Fig. 14) each of the two output parameters, Λ, h, is a function of a single input parameter γ or K: combined with Eqs. (144) and (145). This analytic approximation is compared to the "exact" (numerical) Lagrangian velocity correlation in Fig.  16.

Running Diffusion Coefficient
The running diffusion coefficient for an input algebraic Eulerian correlation is obtained from Eq. (108). Using the algebraic form (111) this integral was calculated analytically in Eqs. (110) and (114). Using also the forms (143) -(145) we obtain the result plotted in Fig. 17 for fixed K = 0.3 and three values of γ = 0.1, 0.3, 0.5. As we know that the regime is always diffusive, the running diffusion coefficient tends toward a finite positive diffusion constant for sufficiently large times (in practice, θ 10). One sees that the asymptotic diffusion coefficient is a weakly decreasing function of the Eulerian exponent.  In most existing works on the Langevin equation the Eulerian velocity time-correlation is assumed to decay exponentially (β = 1). This produces a Lagrangian correlation that also has an exponential asymptotic behavior, and thus automatically ensures the convergence of the integral (108) defining the diffusion constant. It is interesting to compare the diffusion coefficients obtained in the algebraic case and in the exponential case. For the latter we use a Eulerian correlation that is tailored in such a way as to resemble the algebraic one: we thus start the exponential decay at the cut-off time of the algebraic one (Fig. 19): The Lagrangian velocity correlation is calculated in the Corrsin approximation up to the second iteration, by following exactly the same steps as in Secs. 9 and 10. The result is shown in Fig. 20. (Note that L exp (0) = 1, hence it is not necessary to introduce a cut-off time θ L as for the algebraic case). As expected, L exp (θ) decays much more quickly: it has no long-time tail as the algebraic one.
The running diffusion coefficients corresponding to the same parameters are compared in Fig. 21. The algebraic correlation produces a running diffusion coefficient that saturates at a longer time. As a result, the asymptotic diffusion constant D > D exp . This results from the fact that the area under the Lagrangian correlation is larger in the algebraic case (Fig. 20), due to the slow decay of the correlation.

Conclusions
Our main purpose in the present work was to study the relation between a "semi-dynamical" description of a system of particles, based on the V-Langevin equation and its associated hybrid kinetic equation, and on the other hand on a purely stochastic description in terms of a continuous time   random walk and its limiting form as a fractional diffusion equation. We wanted, in particular, to find out under what conditions on the velocity correlation function a V-Langevin equation describes a "strange transport regime", which could be equivalently represented by a CTRW or a FDE.
It was known from previous work [2] that a subdiffusive regime can appear in some limiting cases. This occurs, e.g., for a collisional plasma in a strong magnetic field, in the limit of infinite correlation length perpendicular to B [25], [26] or for electrostatic drift wave turbulence in the limit of infinite Kubo number [2]. In these cases, subdiffusion is produced by sticking of the particles to the magnetic field lines, or by their trapping in the turbulent field. It might be stressed that in the latter case, a diffusive or a subdiffusive (K → ∞) regime was found even when the spatial correlation is long-range (Lorentzian), but the time correlation is exponential [27]. A superdiffusive regime is, however, much more elusive. To the best of our knowledge, a V-Langevin equation producing Lévy-type superdiffusion was only found for a one-dimensional two-state model (velocity = ±V ) [14]. The restriction to a velocity constant in absolute value was necessary in order to replace the Lévy flight CTRW by a Lévy walk.
In the present work we showed that, for one-dimensional systems, the results of [14] concerning the time-fractional superdiffusion can be extended to the more physical case of a general random time-dependent velocity function equipped with an algebraic time-correlation function (but with a Gaussian space-correlation). We confirm that in this case we find a process that is an "intermediate" between diffusion and wave propagation. This type of superdiffusive processes, distinct from the Lévy-type, has been studied in detail in [9] and [10].
We then generalized our investigation to the study of a random velocity field in two dimensions. In this case one is faced with the difficult problem of determining the Lagrangian velocity correlation, a problem that was treated within the Corrsin approximation. It was first shown that if the Eulerian correlation decays algebraically as θ −γ for long times (and its spatial part is Gaussian), the Lagrangian correlation has the same property, with a different exponent: θ −Λ .
We then established a general criterion for the existence of a superdiffusive regime, based on the value of the Lagrangian exponent Λ. It was shown in Sec. 8 that whenever 0 < Λ < 1, the regime is superdiffusive and can be equivalently described by a time-fractional FDE. When Λ > 1, however, the semi-dynamical description predicts a diffusive regime for all Λ, whereas the formally associated FDE predicts subdiffusion, or even looses meaning, for Λ > 2. Thus, the purely stochastic CTRW or FDE descriptions of the turbulence are severely limited in this case by the existence of a bifurcation at Λ = 1.
In the final stage (Sec. 10) the Lagrangian exponent Λ was determined numerically. The result shows that Λ = 2 + γ > 2, for all values of γ and of K. Thus, for the present model and with the approximations used in its treatment: [local approximation of the non-Markovian diffusion equation (54), Corrsin approximation of the Lagrangian velocity correlation (101)]: • The regime is always normal diffusive, for all values of the input Eulerian parameters γ and K. There never appears any superdiffusion, in spite of the long algebraic tail of the temporal velocity correlation. This is to be contrasted with the one-dimensional case of a purely timedependent velocity, where the regime is superdiffusive for all values 0 < γ < 1.
• The "formally associated" fractional differential equation, has no meaning. The evolution cannot be described by a purely stochastic process.
• Although the regime is diffusive, the density profile is not a Gaussian packet. Indeed, although the Lagrangian correlation decays sufficiently fast in order to ensure the existence of a finite asymptotic diffusion constant, the non-Markovian character of the equation (107) cannot be neglected for finite time.
The results obtained here may appear somewhat disappointing, in the sense that no time-fractional superdiffusion was found for two-dimensional motion. Nevertheless, we believe that the result of the exhaustive analysis performed here is interesting in exhibiting the limits of a purely stochastic description. It shows that a long-tailed algebraic temporal correlation is not sufficient for producing superdiffusion (for dimensionality d > 1). In order to obtain a Lévy-type superdiffusive regime, we must presumably act on the spatial part of the Eulerian velocity correlation. This problem is, however, highly non-trivial. On one hand, in the Corrsin approximation, Eq. (106) is invalid, because it is based on the Gaussian form of the spatial Eulerian correlation. But one should go back even further in the chain of approximations. The local approximation leading from Eq. (54) to (55) (in two dimensions) is questionable in the presence of long-range spatial correlations. We hope to come back to this problem in forthcoming work.

Appendix: A Primer on Fractional Calculus
Fractional Calculus has been known to mathematicians for a very long time. The problem appears to have been first raised by L'Hospital and by Leibniz in 1695. The first systematic construction appears in the early 19-th century and is due to Liouville (1832) (well known to every physicist from the first lecture in statistical mechanics), who used this concept in potential theory. It was further developed by B. Riemann in 1847. For a long time the concept was studied only by mathematicians. It entered physics only in the last half of the 20-th century, through the theory of stochastic equations associated with the problems of random walks (Sec. 2).
The search for a natural generalization of the ordinary derivative d n /dx n , with n: an integer, is clearly attractive for a mathematician. It is as natural as the passage from the integers to rational, and later to real numbers. We must find a way to interpolate in a consistent way the concept of differentiation between two successive integers, n and n + 1. The way of doing this is however not obvious. Indeed, several different approaches have been devised by mathematicians. We shall only discuss here the one which is most widely used in the literature (for a discussion of other approaches, see [7], [28] and [11]). One will find detailed expositions of the subject in [8], [13], [7], [29], [5], [6], [30]. Excellent tutorials are the papers by [23] and especially [9], [10]. A brief, but very clear summary is given in [16].
The names of fractional integral, fractional derivative, fractional calculus are actually rather improper. The object of these concepts is to extend the notions of integration and of differentiation from integer order n to arbitrary real order β, not necessarily a fraction. The name has, however, acquired general practice.

RIEMANN-LIOUVILLE FRACTIONAL INTEGRAL
Consider first the multiple integral of order n ( n: integer) 9 of a realvalued function f (x) of the real variable x: The two lower subscripts relate to the limits of integration (a is a fixed real number); the superscript n denotes the order (multiplicity). This equation can be written in the following equivalent form, in terms of a single integral: The equivalence between the two forms is easily demonstrated by partial integration. We may put a J 0 x f (x) = f (x) [7]. For n = 1, the formula reduces to the single integration: Let us check n = 2: In the form (149) the way toward a generalization of the operation of integration (149) is obvious: x a dy f (y) (x − y) 1−α , α > 0 : any positive real number.
(151) This is the definition of the Riemann -Liouville (RL) fractional integral of f (x) of arbitrary real order α.
Some of the more important properties are given here: B) Commutative composition rule: Proof: The ζ-integral is well-known is expressed in terms of the wellknown B-function: hence:

RIEMANN-LIOUVILLE FRACTIONAL DERIVATIVE
The fractional derivative is constructed by combining the operations of the usual (integer-order) derivative with the fractional integral.
We define the fractional derivative as follows: where α is any positive real number, not an integer, and m is the smallest integer larger than α.
In order to "understand intuitively" this rather unusual definition, we may think, very casually, and without any rigour, as follows. We have a solid reference point, viz. a rigorously defined fractional object, the RL integral. We then start with a RL integral of order m − α, where m is the smallest integer which makes m−α > 0, thus ensuring that a J (m−α) x f (x) is a true fractional integral (151). Next, we "differentiate away" the integer m, by taking the m-th ordinary derivative of this object. This looks like "lowering the order" of a J (m−α) x , thus "transforming" it into a J −α x , which could be viewed as a fractional integral of negative order, −α, to be interpreted as a fractional derivative of order +α. But, again, this paragraph is just cavalier talk; in particular, viewing J −α as an "inverse" to D α is an ambiguous statement, as will be shown below. We now go back to business, and consider Eq. (155) as a given definition, and study some of its consequences. A basic requirement is that the definition (155) be compatible with the ordinary (integer-order) derivative. When α is an integer, α = N, we have m − 1 = α = N, thus m = N + 1, and recalling Eq. (150), we have: the operation thus reduces to the ordinary differentiation, as it should. Combining Eqs. (151) and (155) we find: (156) This is the definition of the Riemann -Liouville (RL) left-fractional derivative of order α.
We may also define a Riemann-Liouville (RL) right-fractional derivative in which it is the upper integration limit b that is fixed: We now explore a few properties of the fractional derivatives: we shall meet some surprises. These are all due to a fundamental fact. The fractional derivative (for any of the forms defined above) of a function f (x) is expressed as an integral of f (x) over a given range of x 10 It is therefore an intrinsically non-local operator. In other words, the value of the fractional derivative (for non-integer α) a D α x f (x) at the point x depends on the values of f taken in the whole range of integration defining this quantity. It is quite surprising that, upon varying continuously the order α, this nonlocal character disappears suddenly whenever α reaches any integer value! What is not surprising, therefore, is that fractional calculus entered physics when researchers began to be interested in evolutionary processes endowed with memory ("non-Markovian processes") or in the influence of the spatial neighborhood in the phenomena at a given point (x, t). a) Linearity: The proof follows trivially from the definition of the RL derivative.
b) RL Fractional derivative of a monomial: Proof. Using the definition (156): Let us do this calculation in some detail. Substituting y → xt we find The integral appearing here is expressed in terms of the B-function, (154), hence: Noting the identity: we combine all these partial results to obtain the final formula (159).
The formula (159) is deceptively simple, and would correspond with intuition. Just as the ordinary integral-order derivative (d n /dx n )x m lowers the power of x by n, the fractional derivative transforms x µ into x µ−α , with a prefactor. But a further reflection shows us that this result is not as innocent as it appears. b1) RL Fractional derivative of a constant: Consider, indeed, the case µ = 0, i.e., the i.e., Thus, the RL fractional derivative of non-integer order α > 0 of a constant is non-zero! Is this result consistent with the ordinary derivative? The answer comes from the properties of the gamma function. Indeed, whenever α = n, where n is a non-zero integer, n = 1, 2, 3, ... the argument of the Γ function is a negative integer, or zero. But we know that Hence, for any integer value of α, Eq. (162) reduces, as it should, to: Another feature of Eq. (159) seems to contradict our previous statement: 0 D α x x µ is a perfectly local expression, depending only on the value of x µ at x. 11 In order to understand what happens here, we consider the RL derivative of x µ for a lower integration limit 0 < a < x instead of zero. Starting the calculation as above, we find after the first step: (164) There are two important differences with Eq. (160). First, the integral depends on x, through its lower limit; hence, the m-th derivative acts on the whole product, thus yielding a much more complicated x-dependence. Next, the integral with a non-zero lower limit is no longer a B-function, but rather a very complicated hypergeometric function. In order to illustrate our case, we consider the simplest non-trivial case that leads to simple calculations: µ = 2, α = 1 2 , hence m = 1. Eq. (164) then reduces to: The integral is tabulated; a simple calculation then leads to the result: The result reduces to Eq. (159) when a → 0: The non-local character of the fractional derivative is now luminously illustrated. The complicated function appearing in the right hand side of Eq. (166) is determined by the choice of the function ξ µ in Eq. (164), hence by its values taken over the whole integration range. In particular, the fractional derivative depends explicitly on the lower integration limit a. A striking fact is the singularity of the RL fractional derivative when x → a. This singularity is present for any a, arbitrarily close to zero, but disappears when a = 0. When x ≫ a, a D 1/2 x All this clearly shows that the value a = 0 is not a generic one. The situation is illustrated in Fig. 23. Many authors state rather loosely that "the RL fractional derivative is the inverse of the RL fractional integral". They even adopt the notation: [7]. It will be shown now that this statement is, for the least, ambiguous. Hence Podlubny's notation should preferably be avoided in order to avoid misunderstandings (Mainardi). c1) RL Fractional derivative of a fractional integral of same order.
Proof : For integer α = n ≥ 1: For non-integer α, we enclose this number again between its nearest neighbor integers, defining the integer k as: k − 1 ≤ α < k.
Using now the composition rule (153) for the fractional integrals, we note the identity: Using now the definition (155) we have: (The result follows from the fact that the last step only involves operations of integer order k, i.e., ordinary derivative and integral.) From the result (168) one is tempted to conclude that, indeed, "the RL fractional derivative is the inverse of the RL fractional integral": the fractional RL derivative annuls the action of the fractional integral. The situation is, however, not so clear-cut. It appears from the following property that the double operation D J is not commutative! Indeed: c2) Fractional integral of a RL fractional derivative of the same order.
A partial integration yields: Continuing the partial integrations, and using Eqs. (151) and (153), we end up with: , QED It is precisely because of this non-commutativity of integration and differentiation that one cannot state unambiguously that one operation is the inverse of the other.

d) General composition rules.
In all subsequent formulae we take the lower limit a = 0. We simplify the notations as follows: The proofs of the following results are very similar to the previous ones [7].

A) Fourier transform
The Fourier transform involves integrals ranging from −∞ to +∞. This implies a number of technical assumptions about the regularity of the functions appearing in the theory: we do not discuss them here, assuming that all conditions are satisfied for the existence of the objects of the forthcoming equations. We begin with the following result: Proof: where we set x − y = t. We now note that, from the definition of m, we have: 0 < m − α ≤ 1. We may thus use the tabulated integral: ∞ 0 dξ ξ γ−1 e −ikξ = (ik) −γ Γ(γ), 0 < γ < 1 Hence: This expression is not very convenient in applications. On one hand, the appearance of i α is somewhat awkward. But more important, if we wish to use this result in relation with the Fourier transform, we need to cover the whole range (−∞, ∞), and not only the left range (−∞, x). We therefore consider also the right RL fractional derivative (157): Proof: Setting ξ = −x, we have: Thus: In order to adapt the previous results to the Fourier transform, in particular to cover the whole range (−∞, + ∞), we define a new type of fractional derivative, which is totally symmeric. It is justified by the linearity property of the RL derivative (158): This is called the Riesz fractional derivative. The usefulness of this concept appears obviously in the following result, expressing the Fourier transform of the Riesz derivative: where f (k) denotes the Fourier transform of f (x). Proof: The result (178) ends our quest of the Fourier transform of the (Riesz) fractional derivative of a function, expressed in terms of the Fourier transform of that function.

B) Laplace transform
The (direct) Laplace transform of a function involves an integral over the variable from 0 to ∞. It is therefore particularly well adapted to the study of functions of time, when the latter is restricted to positive values. The two properties of the Laplace transform that are most useful in applications (especially for the solution of differential equations are: The convolution theorem: for any two functions f (t), g(t) of time, defined in the range (0, ∞), their convolution is defined as: Laplace transform of the (ordinary) derivative: where f (n) (t) is the n-th time derivative of f (t), and f (s) is the Laplace transform of f (t). The explicit appearance of the functions f (j) (0) in this formula makes the Laplace transform a particularly useful tool for the solution of initial value problems. We may now note that the value j = 0 in the sum of the right hand side is "special": it corresponds to f (0) (t) = f (t), whereas all other terms in the sum contain true derivatives: We now consider the Laplace transforms of the various objects of fractional calculus.

a) Laplace transform of the fractional integral.
This is easily obtained by noting that the RL integral can be viewed as a convolution: Hence, noting the well-known result: we find immediately, using the convolution theorem (179): b) Laplace transform of the RL derivative.
Note also, from Eq. (184): We thus obtain the final form of the Laplace transform of the RL fractional derivative: .
(191) The analogy with Eq. (181) is striking. The first term has the same form. In the second term the initial value of the function f (0) is replaced, for noninteger β, by the initial value of the fractional integral 0 J −β+m t f (t), whereas in the terms under the sum, the initial values of the ordinary derivatives are replaced by the initial values of the corresponding RL fractional derivatives. In spite of this analogy, the RL time-derivative is not easily adapted to the physical applications. The Laplace transform is well known to be a valuable tool for the solution of initial value problems of differential equations, because it incorporates these initial values in Eq. (181). If we now consider initial value problems of fractional differential equations involving RL fractional derivatives with respect to time, we are faced with a serious difficulty if we try to use a Laplace transformation. The specification of (physical) values of the unknown function and of a sufficient number of its derivatives at t = 0, this does not allow us to determine the correponding initial valued of the fractional integrals and derivatives appearing in Eq. (191). This difficulty has been overcome by M. Caputo in [31], [32] by introducing an alternative definition of the fractional derivative, that is particularly useful for functions of a variable restricted to a finite range, such as time ∈ (0, t).
The Caputo fractional derivative.
The Caputo fractional derivative is defined as follows:  of the RL derivative (191), the Laplace transform of the Caputo derivative only involves the initial values of the function f (t) and of its ordinary derivatives, instead of the initial values of the fractional integral and of the RL fractional derivatives. In any physical differential equation, the former are the given input. This property makes the Caputo derivative an invaluable tool for solving initial value problems of fractional differential equations.