Seismic pulse propagation with constant Q and stable probability distributions

The one-dimensional propagation of seismic waves with constant Q is shown to be governed by an evolution equation of fractional order in time, which interpolates the heat equation and the wave equation. The fundamental solutions for the Cauchy and Signalling problems are expressed in terms of entire functions (of Wright type) in the similarity variable and their behaviours turn out to be intermediate between those for the limiting cases of a perfectly viscous fluid and a perfectly elastic solid. In view of the small dissipation exhibited by the seismic pulses, the nearly elastic limit is considered. Furthermore, the fundamental solutions for the Cauchy and Signalling problems are shown to be related to stable probability distributions with index of stability determined by the order of the fractional time derivative in the evolution equation.

It is known that seismic pulse propagation mostly occurs with a quality factor Q constant over a wide range of frequencies. As pointed out by Caputo and Mainardi (1971) and Caputo (1976), this factor turns out to be independent of frequency only in special linear viscoelastic media for which the stress is proportional to a fractional derivative of the strain, of order ν less than one. Since these media exhibit a creep compliance depending on time by a power-law with exponent ν, we refer to them as power-law solids, according to the notation by Kolsky (1956) and Pipkin (1972Pipkin ( -1986. For the sake of convenience, the generalized operators of integration and differentiation of arbitrary order are recalled in the Appendix in the framework of the so-called Riemann-Liouville Fractional Calculus. In this paper we adopt the Caputo definition for the fractional derivative of order α > 0 of a causal function f (t) (i.e. vanishing for t < 0), where f (m) (t) denotes the derivative of integer order m and Γ is the Gamma function.
In Section 2 we derive the general evolution equation governing the propagation of uniaxial stress waves, in the framework of the dynamical theory of linear viscoelasticity. For a power-law solid the evolution equation is shown to be of fractional order in time, which is intermediate between the heat equation and the wave equation. In fact, denoting the space and time variables by x and t and the response field variable by w(x, t), the evolution equation will be shown to be ∂ 2β w ∂t 2β = D ∂ 2 w ∂x 2 , 2β = 2 − ν . (1. 2) The order of the time derivative has been denoted by 2β for reasons that will appear later. Since 0 < ν ≤ 1 , we get 1/2 ≤ β < 1 .
In Section 3 we review the analysis of the fractional evolution equation (1.2) in the general case 0 < β < 1 , essentially based on our works, Mainardi (1994Mainardi ( , 1995Mainardi ( , 1996aMainardi ( , 1996b. We first analyse the two basic boundary-value problems, referred to as the Cauchy problem and the Signalling problem, by the technique of the Laplace transform and we derive the transformed expressions of the respective fundamental solutions (the Green functions). Then, we carry out the inversion of the relevant Laplace transforms and we outline a reciprocity relation between the Green functions in the space-time domain. In view of this relation the Green functions can be expressed in terms of two interrelated auxiliary functions in the similarity variable r = |x|/( √ Dt β ) . These auxiliary functions can be analytically continued in the whole complex plane as entire functions of Wright type.
In Section 4 we show the evolution of the fundamental solutions for 1/2 ≤ β < 1, that can be relevant in seismology to simulate the propagation of seismic pulses. Accounting for the low dissipation occurring in the Earth, the nearly elastic limit must be considered; in this case the pulse response becomes a narrow, sharply peaked function and the arguments by Pipkin (1972Pipkin ( -1986 and Kreis and Pipkin (1986) must be adopted in order to obtain an evaluation of the solutions, which is suitable from numerical point of view.
Finally, in Section 5, following Kreis and Pipkin (1986), we point out the interesting connection between the fundamental solution for the Signalling problem and the density of a certain (unilateral) stable probability distribution. We note that this connection generalizes the one known for the standard heat equation for which the fundamental solution for the Signalling problem is related to the density of the stable Lévy distribution. Since the above property is expected to provide a further insight into our evolution equation of fractional order, the seismic pulse propagation with constant Q assumes an additional interest from a mathematical-physical point of view.

Linear Viscoelastic Waves and the Fractional Diffusion-Wave Equation
According to the elementary one-dimensional theory of linear viscoelasticity, the medium is assumed to be homogeneous (of density ρ), semi-infinite or infinite in extent (0 ≤ x < +∞ or −∞ < x < +∞) and undisturbed for t < 0 . The basic equations are known to be, see e.g. Hunter (1960), Caputo & Mainardi (1971), Pipkin (1972Pipkin ( -1986, Christensen (1972Christensen ( -1982, Chin (1980), Graffi (1982), Here the suffices x and t denote partial derivation with respect to space and time respectively, the dot ordinary time-derivation, and the star integral timeconvolution from 0 + to t. The following notations have been used: σ for the stress, for the strain, J(t) for the creep compliance (the strain response to a unit step input of stress); the constant J 0 := J(0 + ) ≥ 0 denotes the instantaneous (or glass) compliance.
The evolution equation for the response variable w(x, t) (chosen among the field variables: the displacement u, the stress σ, the strain or the particle velocity v = u t ) can be derived through the application of the Laplace transform to the basic equations. We use the following notation for the Laplace transform of a function f (t) , locally summable for t ≥ 0 , and we adopt the sign ÷ to denote a Laplace transform pair, i.e. f (t) ÷ f (s) .
We first obtain in the transform domain, the second order differential equation is real and positive for s real and positive. As a matter of fact, µ(s) turns out to be an analytic function of s over the entire s-plane cut along the negative real axis; the cut can be limited or unlimited in accordance with the particular viscoelastic model assumed.
Wave like or diffusion like character of the evolution equation can be drawn from (2.5) by taking into account the asymptotic representation of the creep compliance for short times, we have a wave like behaviour with c as the wave-front velocity; otherwise (J 0 = 0) we have a diffusion like behaviour. In the case J 0 > 0 the wave like evolution equation for w(x, t) can be derived by inverting (2.4-5), using (2.6-7) and introducing the non dimensional rate of creep so that the evolution equation turns out to be This is a generalization of D'Alembert wave equation in that it is an integrodifferential equation where the convolution integral can be interpreted as a perturbation term. This case has been investigated by Buchen and Mainardi (1975) and by Mainardi and Turchetti (1975), who have provided wave-front expansions for the solutions.
In the case J 0 = 0 we can re-write (2.6) as where, for the sake of convenience, we have introduced the positive constant D (whose dimensions are L 2 T ν−2 ) and the Gamma function Γ(ν + 1) . Then we can introduce the non-dimensional function φ(t) whose Laplace transform is such that (2.12) Using (2.12), the Laplace inversion of (2.4-5) yields where 2β = 2 − ν so 1/2 ≤ β < 1 . Here the time-derivative turns out to be just the fractional derivative of order 2β (in Caputo's sense), according to the Riemann-Liouville theory of Fractional Calculus recalled in the Appendix.
When the creep compliance satisfies the simple power-law , t > 0 , (2.14) we obtain φ(t) ≡ 0 , and the evolution equation (2.13) simply reduces to (1.2). As pointed out by Caputo and Mainardi (1971), the creep law (2.14) is provided by viscoelastic models whose stress-strain relation (2.3) can be simply expressed by a fractional derivative of order ν . In the present notation this stress-strain relation reads σ = ρD d ν dt ν , 0 < ν ≤ 1 .
(2.15) For ν = 1 the Newton law for a viscous fluid is recovered from (2.15) where D now represents the kinematic viscosity; in this case, since β = 1/2 in (1.2), the classical diffusion equation (or heat equation) holds for w(x, t) . When 0 < ν < 1 the evolution equation (1.2) turns out to be intermediate between the heat equation and the wave equation. In general we refer to (1.2) as the fractional diffusion-wave equation, and its solutions can be interpreted as fractional diffusive waves, see Mainardi (1995).
We point out that the viscoelastic models based on (2.14) or (2.15) with 0 < ν < 1 and henceforth governed by an evolution equation of fractional order in time, see (1.2) with 1/2 < β < 1 , are of great interest in material sciences and seismology. In fact, as shown by Caputo and Mainardi (1971), these models exhibit an internal friction independent on frequency according to the law The independence of the Q from the frequency is in fact experimentally verified in pulse propagation phenomena for many materials including those of seismological interest. From (2.16) we note that the Q is also independent on the material constants ρ and D which, however, play a role in the phenomenon of wave dispersion.
The limiting cases of absence of energy dissipation (the elastic energy is fully stored) and of absence of energy storage (the elastic energy is fully dissipated) are recovered from (2.16) for ν = 0 (perfectly elastic solid) and ν = 1 (perfectly viscous fluid), respectively.
To obtain values of seismological interest for the dissipation (Q ≈ 1000) we need to choose the parameter ν sufficiently close to zero, which corresponds to a nearly elastic material; from (2.16) we obtain the approximate relations between ν and Q , namely As a matter of fact the evolution equation (1.2) turns out to be a linear Volterra integro-differential equation of convolution type with a weakly singular kernel of Abel type. Equations of this kind have been treated, both with and without reference to the fractional calculus, by a number of authors including Caputo (1969Caputo ( , 1976Caputo ( , 1996b, Meshkov and Rossikhin (1970), Pipkin (1972Pipkin ( -1986, Buchen and Mainardi (1975), Kreis and Pipkin (1986), Nigmatullin (1986), Schneider and Wyss (1989), Giona and Roman (1992), Metzler et al. (1994) and Mainardi (1994Mainardi ( , 1995Mainardi ( , 1996aMainardi ( , 1996b. For recent reviews on related topics we refer to Rossikhin and Shitikova (1997) and Mainardi (1997).

The Reciprocity Relation and the Auxiliary Functions
The two basic problems for our fractional wave equation (1.2) concern, for t ≥ 0, the infinite interval −∞ < x < +∞ and the semi-infinite interval x ≥ 0 , respectively; the former is an initial -value problem, referred to as the Cauchy problem, the latter is an initial boundary -value problem, referred to as the Signalling problem.
Extending the classical analysis to our fractional equation (1.2), and denoting by g(x) and h(t) two given, sufficiently well-behaving functions, the basic problems are thus formulated as following, a) Cauchy problem, If 1/2 < β < 1 , we must add in (3.1a) and (3.1b) the initial values of the first time derivative of the field variable, w t (x, 0 + ) , since in this case (1.2) contains a time derivative of the second order. To ensure the continuous dependence of our solution with respect to the parameter β also in the transition from β = (1/2) − to β = (1/2) + , we agree to assume w t (x, 0 + ) = 0 .
In view of our analysis we find it convenient from now on to add the parameter β to the independent space-time variables x , t in the solutions, writing w = w(x, t; β) .
For the Cauchy and Signalling problems we introduce the so-called Green functions G c (x, t; β) and G s (x, t; β), which represent the respective fundamental solutions, obtained when g(x) = δ(x) and h(t) = δ(t) . As a consequence, the solutions of the two basic problems are obtained by a space or time convolution according to since the Green function of the Cauchy problem turns out to be an even function of x. According to a usual convention, in (3.2b) the limits of integration are extended to take into account for the possibility of impulse functions centred at the extremes.
For the standard diffusion equation (β = 1/2) it is well known that In the limiting case β = 1 we recover the standard wave equation, for which, In the general case 0 < β < 1 the two Green functions will be determined by using the technique of the Laplace transform. This technique allows us to obtain the transformed functions G c (x, s; β), G s (x, s; β), by solving ordinary differential equations of the 2-nd order in x and then, by inversion, G c (x, t; β) and G s (x, t; β).
For the Cauchy problem (3.1a) the application of the Laplace transform to (1.2) with w(x, t) = G c (x, t; β) leads to the non homogeneous differential equation satisfied by the image of the Green function, G c (x, s; β) , Because of the singular term δ(x) we have to consider the above equation separately in the two intervals x < 0 and x > 0, imposing the boundary conditions at x = ∓∞ , G c (∓∞, t; β) = 0 , and the necessary matching conditions at x = 0 ± . We obtain For the Signalling problem (3.1b) the application of the Laplace transform to (1.2) with w(x, t) = G s (x, t; β) leads to the homogeneous differential equation Imposing the boundary conditions at From (3.6) and (3.8) we recognize for the original Green functions the following reciprocity relation This relation can be easily verified in the case of standard diffusion (β = 1/2), where the explicit expressions (3.4) of the Green functions leads to the identity We refer to F d (r) and M d (r) as to the auxiliary functions for the diffusion equation because each of them provides the fundamental solution through (3.10). We note that M d (r) satisfies the normalization condition ∞ 0 M d (r) dr = 1. Applying in the reciprocity relation (3.9) the complex inversion formula for the transformed Green functions (3.6) and (3.8), and changing the integration variable in σ = s t , we obtain ( 3.12) where is the similarity variable and F (r; β) := 1 2πi Br e σ − rσ β dσ , M (r; β) := 1 2πi Br e σ − rσ β dσ σ 1−β (3.14) are the auxiliary functions. In (3.14) Br denotes the Bromwich path and r > 0 , 0 < β < 1 .
The above definitions of F (r; β) and M (r; β) by the Bromwich representation can be analytically continued from r > 0 to any z ∈ C, by deforming the Bromwich path Br into the Hankel path Ha , a contour that begins at σ = −∞ − ia (a > 0), encircles the branch cut that lies along the negative real axis, and ends up at σ = −∞ + ib (b > 0).
The integral and series representations of F (z; β) and M (z; β), valid on all of C , with 0 < β < 1 turn out to be In the theory of special functions, see Erdélyi (1955), we find an entire function, referred to as the Wright function, which reads (in our notation) Although convergent in all of C, the series representations in (3.15-16) can be used to provide a numerical evaluation of our auxiliary functions only for relatively small values of r , so that asymptotic evaluations as r → +∞ are required. Choosing as a variable r/β rather than r , the computation by the saddle-point method for the M function is easier and yields, see Mainardi and Tomirotti (1995), We note that the saddle-point method for β = 1/2 provides the exact result (3.11), i.e. M (r; 1/2) = M d (r) = exp(−r 2 /4)/ √ π , but breaks down for β → 1 − . The case β = 1 , for which (1.2) reduces to the standard wave equation, is of course a singular limit also for the series representation since M (r; 1) = δ(r − 1).

The Evolution of the Seismic Pulse from the Fundamental Solutions
It is known that in theoretical seismology the delta-Dirac function is of great relevance in simulating the pulse generated by an ideal seismic source, concentrated in space (δ(x)) or in time (δ(t)). Consequently, the fundamental solutions of the Cauchy and Signalling problems are those of greater interest because they provide us with information on the possible evolution of the seismic pulses during their propagation from the seismic source. Accounting of the reciprocity relation (3.12) and the similarity variable (3.13), r = x/( √ D t β ) , the two fundamental solutions can be written, for x > 0 and t > 0 , as (4.1b) The above equations mean that for the fundamental solution of the Cauchy [Signalling] problem the time [spatial] shape is the same at each position [instant], the only changes being due to space [time] -dependent changes of width and amplitude. The maximum amplitude in time [space] varies precisely as 1/x [1/t]. The two fundamental solutions exhibit scaling properties that make easier their plots versus distance (at fixed instant) and versus time (at fixed position). In fact, using the well-known scaling properties of the Laplace transform in (3.6) and (3.8), we easily prove, for any p , q > 0 , that and, consequently, in plotting we can choose suitable values for the fixed variable. In order to inspect the evolution of the initial pulse for seismological purposes, we need to plot G c (x, t; β) versus x and G s (x, t; β) versus t as β is sufficiently close to 1 (nearly elastic cases) to ensure a sufficiently low value for the constant internal friction Q −1 . From (2.13)-(2.14) and (2.16)-(2.17) we need to consider β = 1− with = ν/2 of the order of 0.001 to 0.01 . In the evaluation of the auxiliary functions in the nearly elastic cases, we note that the matching between the series and saddle point representations is no longer achieved since the saddle point turns out to be wide and the consequent approximation becomes poor. In these cases we need to adopt the ingenious variant of the saddle-point method introduced by Pipkin (1972Pipkin ( -1986, see also Kreiss and Pipkin (1986), which allows us to see some structure in the peak while it tends to the Dirac delta function. With Pipkin's method we get the desired matching with the series representation just in the region around the maximum r ≈ 1 , as shown in Fig. 2a,b, where we exhibit the significant plots of the auxiliary function M (r; β) with β = 1 − for = 0.01 and = 0.001 .  We also note the exponential decay of G c (x, t; β) as x → +∞ (at fixed t) and the algebraic decay of G s (x, t; β) as t → +∞ (at fixed x), for 0 < β < 1 . In fact, using (4.1a,b) with (3.16) and (3.19), we get where a(t) , b(t) and c(x) are positive functions.

The Fundamental Solutions as Probability Density Functions
It is well known that the fundamental solution of the standard diffusion equation for the Cauchy problem is related with the Gauss or normal probability law, bilateral in space. In fact, recalling (3.3a), we have denotes the well-known Gauss probability density function (pdf ) with variance σ 2 . The associated cumulative distribution function (cdf ) is known to be The moments of even order of the Gauss pdf turn out to be with n ∈ IN, If we consider the fundamental solution of the standard diffusion equation but for the Signalling problem, we note that it is related to the Lévy probability law, unilateral in time (a property not so well-known as that for the Cauchy problem!). In fact, recalling (3.3b), we have denotes the Lévy pdf , see Feller (1971), with cdf The Lévy pdf has all moments of integer order infinite, since it decays at infinity as t −3/2 . However, we note that the moments of real order δ are finite only if 0 ≤ δ < 1/2 . In particular, for this pdf the mean (expectation) is infinite, but the médiane is finite. In fact, from P L (t med ; µ) = 1/2 , it turns out that t med ≈ 2µ , since the complementary error function gets the value 1/2 as its argument is approximatively 1/2 (a better evaluation of the argument is 1/2.1).
The Gauss and Lévy laws are special cases of the important class of αstable probability distributions, or stable distributions with index of stability (or characteristic exponent) α = 2 and α = 1/2 , respectively. Another special case is provided by the Cauchy law with pdf p C (x; λ) = λ/[π(x 2 + λ 2 )] and α = 1 .
The name stable has been assigned to these distributions because of the following property: if two independent real random variables with the same shape or type of distribution are combined linearly and the distribution of the resulting random variable has the same shape, the common distribution (or its type, more precisely) is said to be stable. More precisely, if Y 1 and Y 2 are random variables having such distribution, then Y defined by the linear combination c Y = c 1 Y 1 + c 2 Y 2 has a similar distribution with the same index α for any positive real values of the constants c , c 1 and c 2 with c α = c α 1 + c α 2 . As a matter of fact only the range 0 < α ≤ 2 is allowed for the index of stability. The case α = 2 is noteworthy since it corresponds to the normal distribution, which is the only stable distribution which has finite variance, indeed finite moments of any order. In the cases 0 < α < 2 the corresponding pdf p α (y) have inverse power tails, i.e. |y|>λ p α (y) dy = O(λ −α ) and therefore their absolute moments of order δ are finite if 0 ≤ δ < α and infinite if δ ≥ α .
The inspiration for systematic research on stable distributions, originated with Paul Lévy, was the desire to generalize the celebrated Central Limit Theorem (CLT ).
The restrictive condition of stability enabled some authors to derive the general form for the characteristic function (cf , the Fourier transform of the pdf ) of a stable distribution, see Feller (1971). A stable cf is also infinitely divisible, i.e. for every positive integer n it can be expressed as the nth power of some cf . Equivalently we can say that for every positive integer n a stable pdf can be expressed as the n-fold convolution of some pdf . All stable pdf are unimodal and indeed bell-shaped, i.e. their n-th derivative has exactly n zeros, Using standardized random variables, the α-stable distributions turn out to depend on an additional parameter γ , said the skewness parameter. Denoting a stable pdf by p α (y; θ) , we note the property p α (−y; −θ) = p α (y; θ) . Consequently a stable pdf with θ = 0 is necessarily symmetrical. As a matter of fact |θ| ≤ α if 0 < α < 1 and |θ| ≤ 2 − α if 1 < α < 2 . Stable distributions with extremal values of θ are called extremal.
We easily recognize that for 0 < β < 1 the fundamental solution for the Signalling problem provides a unilateral extremal stable pdf in (scaled) time with index of stability α = β , which decays according to (4.3b) with a power law. In fact, from (4.1b) and (5.11) we note that, putting y = r −1/β = τ , This property has been noted also by Kreiss and Pipkin (1986) based on (3.8) and on Feller's result, p α (t; −α) ÷ exp(−s α ) for 0 < α < 1 .
As far as the Cauchy problem is concerned, we note that the corresponding fundamental solution provides a bilateral symmetrical pdf in (scaled) distance with two branches, for x > 0 and x < 0 , obtained one from the other by reflection. For large |x| each branch exhibits an exponential decay according to (4.3) and, only for 1/2 ≤ β < 1 , it is the corresponding branch of an extremal stable pdf with index of stability α = 1/β . In fact, from (4.1b) and (5.12) we note that, putting y = |r| = ξ > 0 , (5.14) This property had to the authors' knowledge not been noted: it properly generalizes the Gaussian property of the pdf found for β = 1/2 (standard diffusion). Furthermore,using (3.20), the moments (of even order) of G c (x, t; β) turn out to be We recognize that the variance is now proportional to Dt 2β , which implies a phenomenon of fast diffusion if 1/2 < β < 1 .

Appendix: Essentials of Fractional Calculus
Fractional calculus is the field of mathematical analysis which deals with the investigation and applications of integrals and derivatives of arbitrary order. The term fractional is a misnomer, but it is retained following the prevailing use.
According to the Riemann-Liouville approach to fractional calculus, the notion of fractional integral of order α (α > 0) is a natural consequence of the well known formula (usually attributed to Cauchy), that reduces the calculation of the n−fold primitive of a function f (t) to a single integral of convolution type. In our notation the Cauchy formula reads where IN is the set of positive integers. From this definition we note that f n (t) vanishes at t = 0 with its derivatives of order 1, 2, . . . , n − 1 . For convention we require that f (t) and henceforth f n (t) be a causal function, i.e. identically vanishing for t < 0 .
In a natural way one is led to extend the above formula from positive integer values of the index to any positive real values by using the Gamma function. Indeed, noting that (n − 1)! = Γ(n) , and introducing the arbitrary positive real number α , one defines the Fractional Integral of order α > 0 : where IR + is the set of positive real numbers. For complementation we define J 0 := I (Identity operator), i.e. we mean J 0 f (t) = f (t) . Furthermore, by J α f (0 + ) we mean the limit (if it exists) of J α f (t) for t → 0 + ; this limit may be infinite. We note the semigroup property J α J β = J α+β , α , β ≥ 0 , which implies the commutative property J β J α = J α J β , and the effect of our operators J α on the power functions These properties are of course a natural generalization of those known when the order is a positive integer. Introducing the Laplace transform by the notation L {f (t)} := ∞ 0 e −st f (t) dt = f (s) , s ∈ C , and using the sign ÷ to denote a Laplace transform pair, i.e. f (t) ÷ f (s) , we note the following rule for the Laplace transform of the fractional integral, which is the generalization of the case with an n-fold repeated integral. After the notion of fractional integral, that of fractional derivative of order α (α > 0) becomes a natural requirement and one is attempted to substitute α with −α in the above formulas. However, this generalization needs some care in order to guarantee the convergence of the integrals and preserve the well known properties of the ordinary derivative of integer order.
Denoting by D n with n ∈ IN , the operator of the derivative of order n , we first note that D n J n = I , J n D n = I , n ∈ IN , i.e. D n is left-inverse (and not right-inverse) to the corresponding integral operator J n . In fact we easily recognize from (A.1) that As a consequence we expect that D α is defined as left-inverse to J α . For this purpose, introducing the positive integer m such that m − 1 < α ≤ m , one defines the Fractional Derivative of order α > 0 as Defining for complementation D 0 = J 0 = I , then we easily recognize that D α J α = I , α ≥ 0 , and Of course, these properties are a natural generalization of those known when the order is a positive integer. Note the remarkable fact that the fractional derivative D α f is not zero for the constant function f (t) ≡ 1 if α ∈ IN . In fact, (A.7) with γ = 0 teaches us that This, of course, is ≡ 0 for α ∈ IN, due to the poles of the gamma function in the points 0, −1, −2, . . .. We now observe that an alternative definition of fractional derivative, originally introduced by Caputo (1967Caputo ( ) (1969 in the late sixties and adopted by Caputo and Mainardi (1971) in the framework of the theory of Linear This definition is of course more restrictive than (A.6), in that requires the absolute integrability of the derivative of order m. Whenever we use the operator D α * we (tacitly) assume that this condition is met. We easily recognize that in general unless the function f (t) along with its first m − 1 derivatives vanishes at t = 0 + . In fact, assuming that the passage of the m-derivative under the integral is legitimate, one recognizes that, for m − 1 < α < m and t > 0 , .11) and therefore, recalling the fractional derivative of the power functions (A.7) , The alternative definition (A.9) for the fractional derivative thus incorporates the initial values of the function and of its integer derivatives of lower order. The subtraction of the Taylor polynomial of degree m − 1 at t = 0 + from f (t) means a sort of regularization of the fractional derivative. In particular, according to this definition, the relevant property for which the fractional derivative of a constant is still zero can be easily recognized, i.e.
We now explore the most relevant differences between the two fractional derivatives (A.6) and (A.9). We agree to denote (A.9) as the Caputo fractional derivative to distinguish it from the standard Riemann-Liouville fractional derivative (A.6). We observe, again by looking at (A.7), that D α t α−1 ≡ 0 , α > 0 , t > 0 . From above we thus recognize the following statements about functions which for t > 0 admit the same fractional derivative of order α , with m − 1 < α ≤ m , m ∈ IN , D α f (t) = D α g(t) ⇐⇒ f (t) = g(t) + m j=1 c j t α−j , (A.14) D α * f (t) = D α * g(t) ⇐⇒ f (t) = g(t) + In these formulas the coefficients c j are arbitrary constants.
For the two definitions we also note a difference with respect to the formal limit as α → (m − 1) + . From (A.6) and (A.9) we obtain respectively, We now consider the Laplace transform of the two fractional derivatives. For the standard fractional derivative D α the Laplace transform, assumed to exist, requires the knowledge of the (bounded) initial values of the fractional integral J m−α and of its integer derivatives of order k = 1, 2, . . . , m − 1 . The corresponding rule reads, in our notation, The Caputo fractional derivative appears more suitable to be treated by the Laplace transform technique in that it requires the knowledge of the (bounded) initial values of the function and of its integer derivatives of order k = 1, 2, . . . , m − 1 , in analogy with the case when α = m . In fact, by using (A.4) and noting that we easily prove the following rule for the Laplace transform, Indeed, the result (A.20), first stated by Caputo (1969) by using the Fubini-Tonelli theorem, appears as the most "natural" generalization of the corresponding result well known for α = m . This appendix is based on the review by Gorenflo and Mainardi (1997). For more details on the classical treatment of fractional calculus the reader is referred to Erdélyi (1954), Oldham and Spanier (1974), Samko et al. (1987Samko et al. ( -1993 and Miller and Ross (1993). Gorenflo and Mainardi have pointed out the major utility of the Caputo fractional derivative in the treatment of differential equations of fractional order for physical applications. In fact, in physical problems, the initial conditions are usually expressed in terms of a given number of bounded values assumed by the field variable and its derivatives of integer order, no matter if the governing evolution equation may be a generic integro-differential equation and therefore, in particular, a fractional differential equation.