About Some Possible Implementations of the Fractional Calculus

: We present a partial panoramic view of possible contexts and applications of the fractional calculus. In this context, we show some different applications of fractional calculus to different models in ordinary differential equation (ODE) and partial differential equation (PDE) formulations ranging from the basic equations of mechanics to diffusion and Dirac equations.


From Elementary Mathematical Analysis to Fractional Derivatives
A first reason to justify the use of fractional operators is the need to introduce memory terms into differential models in a natural form. In this sense, we can consider the classical Calculus Fundamental Theorem with the introduction of a convolution kernel F associated to a function g leads to the natural form of introducing a memory term by changing the convolution kernel of the integral: F(x) = F(a) + x a K(g(x) − g(s)) · f (s)ds (1) Now F is the generalized primitive of the function f and the convolution kernel K is a memory term that could be different in each specific problem, changing the definition of fractional operator consequently. This generalization of the integral can be considered as a base to construct possible definitions for fractional integrals.
From these considerations, the fractional calculus emerges in the mathematical world as the study of integral and derivative operators of non-integer orders on domains of real or complex functions. Several definitions of fractional derivatives D α have been developed progressively with the objective to generalize the concept of ordinary derivative D, such that for α = 1 the ordinary operator can be recovered [1][2][3].
Fractional calculus is a powerful mathematical tool that allows to create intermediate-order parameters equations and offers modeling scenarios where fundamental mathematical questions converge and appropriate numerical algorithms can de developed. A lot of fractional operators have been defined in the literature; however, not all of them can be used in each real-world application. In this context, we appreciate very much the enthusiastic and clarifier paper of D. Baleanu and A. Fernandez [4] and M.D. Ortigueira and J.A.T. Machado [5][6][7].

From Factorial to the Gamma Function
The first definitions of fractional operators are related to the use of the Gamma function is a function that generalizes the definition of the factorial to non-positive numbers. Its definition is: for any complex number z with positive real part. Using integration by parts in Equation (2), a fundamental property of the Gamma function is obtained: which allows to give the Gamma function of a positive integer number as Γ(n) = (n − 1)! . (4) In this context, the Gamma function is a generalization of the concept of factorial.
In 1738, Euler introduces the first generalization of ordinary derivative, verifying that the fractional derivation made sense for the potential function x a . And in 1819, Lacroix starts from the m-order derivative of the function x n , with m and n positive integer numbers to determine the 1/2 order derivative of the function x a , using the generalization of the factorial function by the Gamma function: such that for a = 1: This is the result that will be obtained with the called Riemann-Lioville fractional derivative. This concept can be generalized to any order and the following relation between the ordinary and fractional case with the Riemann-Liouville fractional derivative is obtained: Later, some provisional definitions of fractional operators were introduced by Fourier, Abel, Liouville and Riemman, without much success. Until, in 1870, N. Ya. Sonine started from the Cauchy formula for repeated integration: x 1 a dx 2 · · · x n−1 a f (t)dt = 1 (n − 1)! x a (x − t) n−1 f (t)dt (9) and, using the generalization of the factorial function by the Gamma function, he obtained the actual definition of the fractional integral of Riemann-Liouville: although in 1884 Laurent formulated it definitively.

Some Definitions of Fractional Integrals and Derivatives
An important definition of fractional integral and derivative corresponds to Riemann-Liouville: • Left-side Riemann-Liouville Fractional Integral of order α > 0: • Right-side Riemann-Liouville Fractional Integral of order α < 0: • Left-side Riemann-Liouville Fractional Derivative of order α > 0: • Right-side Riemann-Liouville Fractional Derivative of order α < 0: In all cases n ∈ N, such that 0 ≤ n − 1 < α < n. These operators recover the classical operators for the parameter α = 1 and the algebra of these operators is different to the classical operators: and Re(α), Re(β) > 0. Then: . Then: Related to the integrals of Riemann-Liouville, the definition of Caputo fractional derivative appears: • Left-side Caputo Fractional Derivative of order α > 0: • Right-side Caputo Fractional Derivative of order α < 0: We have, as before, n ∈ N such that 0 ≤ n − 1 < α < n, and now the n + 1 derivatives of function φ must be continuous and bounded in [a, b].
The following identity established the relation between the Riemann-Liouville and Caputo fractional derivatives, for f a suitable function (for instance, f n-derivable): The extension of Riemann-Liouville fractional operators to infinity intervals leads to the Liouville fractional operators: • Left-side Liouville Fractional Integral of order α > 0: • Right-side Liouville Fractional Integral of order α < 0: • Left-side Liouville Fractional Derivative of order α > 0: • Right-side Liouville Fractional Derivative of order α < 0: (25) In all cases n ∈ N, such that 0 ≤ n − 1 < α < n.

Mittag-Leffler Functions
Special functions related to the eigenfunctions of fractional operators are the Mittag-Leffler functions. They appear in the solution of many fractional differential equations. The Mittag-Leffler functions are generalizations of the exponential function and they was introduced by the mathematician G.M. Mittag-Leffler in 1903: Here, we have some elementary properties of the Mittag-Leffer function. For some values of the parameters α, β, Mittag-Leffer functions return known classical functions, for example: The relevance of the Mittag-Lefler functions is their behavior as generalized exponential functions associated to the Riemann-Liouville and Caputo fractional derivatives: On the other hand, a very useful extension of the Mittag-Leffer function is related to extend the fractional derivatives defined by using a Mittag-Leffler kernel which is non-local and non-singular ( [8,9]).
valid for 0 < α < 1, with B(α) being a normalisation function. This new definition has many applications at the same time that satisfies the extensions of the product rule and chain rule. Futhermore, these functions have many uses, for instance, they allow to address fractal kinetics from basic functions ( [10]) or they can be used as a simplification tool combined to exponent/powerlaw mathematical congruence ( [11]).

Some Ideas for Numerical Integration
It is simple to extend some classical methods from integer to non-integer orders. For instance, the most basic first order explicit method for numerical integration of ordinary differential equations with a given initial value, known as Euler methods.
In classical dynamical models, symplectic integration schemes preserve the flow of the hamiltonian while other classical integrators as Runge-Kutta schemes do not necessarily conserve it. This is immediately translated into the conservation of the first integral of motion and long term stability of the scheme. In fractional dynamics, energy is not conserved. Despite this, using a symplectic scheme in fractional mechanics ensures that the observed instabilities will be certainly due to the fractional operators. Long-term stability is inherited in this fractional mapping as is shown below.
Let us consider the following initial value problem with Caputo fractional derivative: where y α is a Lipschitz function with constant L. By using a truncated Taylor series for the Caputo fractional derivative, with step t n = nh with h = T N y n = 0, 1, 2, . . . , N, we obtain: and then the fractional Euler method for this initial value problem is: f (t n , y n ) (36) with convergence order h α . Other more complex generalizations of classical methods are possible. For instance, fractional Hamiltonian-Jacobi methods. The equations of motions for a one-dimensional Hamiltonian system H = 1 2 p 2 + V(x) with unit mass, are defined aṡ that is associated to the second order equationẍ + V (x) = 0. This system can be generalized by using Caputo time-fractional derivatives where 0 < α ≤ 1.
And the system Equation (38) with initial condition x(0) = x 0 , p(0) = p 0 is equivalent to For the numeric solution of system Equation (39) we have developed a map (see [12]) When α = 1, this is equivalent to a second order symplectic integrator p n = p n−1 − ∆t V (x n−1 ), x n = x n−1 + ∆t p n , and mapping Equation (38) provides an orbit (x n , p n ) approaching the exact orbit at t = n∆t when ∆t → 0. Futhermore, the term p k+1 can be replaced by p k in order to return the second order Euler scheme which is not symplectic as ∆t → 0.
The orbit at step n depends on all the previous states up to the initial one due to the memory kernel of the fractional integral and then we have an infinite dimensional mapping. The computational complexity of the orbit up to (x n , p n ) is of order n 2 whereas it is of order n for α = 1.
This map has been tested taking α = 1 with standard models and using different potentials which solutions are known. In particular for the harmonic oscillator and initial conditions x 0 = 1, p 0 = 0 and ∆t = 0.01, solution has been compared with x(t) = cos(t), p(t) = sin(t) with error smaller than 0.5% after 10,000 steps.
The map Equation (40) has been used to simulate numerical solutions to significant non-linear fractional generalized Hamiltonian problems with a potential V(x), where the explicit solution cannot be found. For instance, standard academic cases like free particle motion (V = 0) and a uniformly accelerated particle (V = kx) [13], or the simple oscillator V = 1 2 ω 2 x 2 , the double well potential Figure 1) and the pendulum V = cos(x) [12]. In this pioneer numerical work, it was observed that the fractional derivative introduces a damping effect which can be either algebraic or exponential depending on the time scales of the system. It is considered in a more general context in [10] or in [11].
Riemann-Liouville time fractional problem could be integrated in this way, changing Riemann-Liouville derivative by Caputo's and applying a similar mapping. Non-homogeneous systems can be studied through the generalizations of this map. For instance, by the introduction of an external force f (t) to simulate a forced-damped oscillator [14]: and this change means the introduction of an extra force term f n = f (n∆t) in the first equation of (40).
In particular, the system evolves to a limit cycle for an harmonic forcing ( f (t) = A 0 cos(ωt)), similar to the classic case with the forced-damped oscillator. Varying the forcing frequency ω a resonance motion is reproduced, with amplitude (see Figure 2).

Nonlocal, Fractional Calculus of Variations
Classical mechanics can be viewed as founded on Hamilton's principle of stationary action. It is an elegant theory based on a very simple axiom. Some authors, appealed both by this and by the potential effectiveness of fractional modeling, have developed a fractional mechanics.
The guiding idea is to keep the same axiom (Hamilton's principle) and allow to have fractional derivatives as variables. Apparently, the founding papers are due to Riewe [15,16]. He uses fractional derivatives and builds the corresponding Euler-Lagrange equations in a systematic way. He also gives the corresponding Hamilton equations. As an application, he provides a formulation that includes a dissipative force. His presentation is not without some limitation.
Agrawall and other authors have generalized from them a Lagrangian and a Hamiltonian formalism [17][18][19] that also includes constraints.
From these references we see that right and left fractional derivatives appear in the Lagrangian and in the fractional Euler-Lagrange equations. The implications that this poses to the causality have not been dealt with.
We propose a slightly more general formalism that allows systems with different nonlocal terms in the Lagrangian, including fractional integrals and fractional derivatives of Caputo and Riemann-Liouville kind.

Positions
Let us call t, the time, the independent variable, x(t) a function and y some linear functional that depends on x over the whole range of times through a relation yet to be defined: We consider a (Lagrangian) function of these, L(t, x, y), and we look for the corresponding E-L equation for the action L using the ideas of the calculus of variation. For instance, following Agrawall [17], we consider x * (t) and y * (t) the functions that make the action stationary and write x(t) = x * (t) + εη(t). This implies, due to the assumed linearity of G that With this, the action becomes a function of ε with derivative: The next step is to write the integrand in terms of η(t) and declare each factor to be zero. For this, we need to suppose a specific form to the functional G.
Let us suppose a nonlocal dependence (in time) of y on x of the form where K(t, s) is some kernel, independent of both x and y. With this we have from Equation (45) T 0 The use of Dirac's delta is justified if we understand that L is zero for t outside [0, T]. Besides, for all this manipulation to make sense, we need the kernel K and both partials of L to be integrable and ∂L/∂x to have a finite L 2 -norm. If we choose a singular kernel, as below, we will see that some other conditions might be necessary.
Equation (47) is the corresponding E-L equation, we may exchange in it the name of the variables t and s and we finally have as necessary condition for the stationary action: If we compare the integral with Equation (46) we see that the variables in the kernel are interchanged.
This is an important feature and it implies that the equation at a given time t involves values of x from both the past and the future. We will see this more clearly below when we deal with fractional derivatives.

Velocities
Let us suppose now that, instead of depending on x, y(t) depends linearly onẋ(s) through We repeat what we did previously: we consider x * (t), y * (t) and write With this, the action becomes again a function of ε with the condition of the stationary trajectories given by: Substituting the value of G, we have From here we get as Euler-Lagrange equation where we have interchanged in the last step the names of the variables s and t. As in the previous case, where y depended only on positions, the evolution of the system at time t depends on both past and future values.

Extension of the Velocities Case
We may also consider the following possibility: y depends on x as given by Equation (46), and L depends on x and y but also onẏ.
Since the procedure is linear, we may consider just that L depends on x andẏ. The necessary condition for stationary trajectories becomes: Operating similarly as before, we have: Or, otherwise: if we consider that L and its derivatives are all zero at the boundaries. As before, we may rewrite this as: or, as: We have, once more, this dependency on values from the past and from the future.

General Case
We may consider a Lagrangian with three variables, for instance, y 1 depending on x,ẏ 1 , and y 2 depending onẋ, through two linear functionals, as above, with kernels K 1 and K 2 , respectively: In that case, due to the linearity of all the previous manipulations, the Euler-Lagrange equation will just have a contribution from each and be of the form: This reminds us of the original calculus of variations (i.e., local and nonfractional) with higher order formulation. The case, for instance, of a Lagrangian that depends on x,ẋ andẍ.

Fractional Integrals
We may obtain fractional integrals for y given by Equation (46) choosing, for instance, the following kernel: where α > 0 and Θ is the heavyside function: Kernel Equation (61) gives us the left-sided Riemann-Liouville integral with lower boundary 0: while interchanging the variables and considering K I − (t, s) = K I + (s, t), we obtain the right-sided integral with upper boundary T: For y(t) = 0 + I α t x(t), the Euler-Lagrange Equation (48) becomes: while in the second case, where y(t) = T − I α t x(t), we obtain: We see that, independently of the choice we consider, the Euler-Lagrange equation for the system has both kinds of integrals: left-sided and the right-sided, one as the variable y, the other applied to the partial derivative of L with respect to y. This supposes that the evolution equation at any intermediate time t ∈ (0, T), has elements from both the past and the future.

Fractional Derivatives
Caputo: In order to have y to represent a fractional derivative, we may consider the dependence on velocities and use in Equation (49) the kernel We have that y(t) is the (left) fractional Caputo derivative of order α: The Euler-Lagrange equation is in this case Conversely, if we choose in Equation (49) the kernel we have that y is the right Caputo derivative of order α and the corresponding Euler-Lagrange Equation (53) is now We see that for both kernels we have in the Euler-Lagrange equations derivatives form both sides, as in the previous case.
In both circumstances we have an equation that involves values from the past and from the future.
if we take the left-derivative for y, and is: if we chose the right-derivative. By the way, the sign that seemed to be missing ("the force is minus the derivative of the potential") is included inside the right-derivative in both cases, as we have seen right above for one of them. (46) the kernel:

Riemann-Liouville: we have to use the extension of the velocities case and consider in Equation
since it yields:ẏ The corresponding Euler-Lagrange equation is: Conversely, if we consider the kernel: we obtain:ẏ and the Euler-Lagrange equation: As always, we obtain integrals that cover both ranges: from 0 to t and from t to T.
We also see that the result is very similar that for the Caputo derivative, with an exchange of the role of the Caputo and the Riemann-Liouville derivatives.

Momenta and Hamilton Formalism
For each variable y we may define a momentum associated to it, in analogy to classical mechanics, as or and we can express the previous Euler-Lagrange equations as ∂L(t, x, y) ∂x when y depends on positions, and when depending on velocities (first case) or (second case) Using the four kernels considered before and the corresponding Euler-Lagrange equations we may give the Hamiltonian approach to the systems.
For instance, when y depends on velocities, the momentum is given by Equation (80), we have with kernel Equation (68): and with kernel Equation (71): When y depends on positions, the Lagrangian depends onẏ and the momentum is Equation (81), we obtain with kernel Equation (74): and with kernel (77): Whenever from Equation (80) we can express y (orẏ) as a function of x and p, i.e., we can provide a Legendre transformation, the previous systems of equations correspond to the Hamilton equations.
As for the general case presented in Section 2.1.4, we may consider several functions y in the Lagrangian, and each will give rise to its corresponding momentum and a corresponding evolution equation for it.

Interpretation
We have seen that nonlocality (as considered above) plus Hamilton's principle implies that the equations are not causal, in the sense that the evolution at a given time t depends on values from all times and not just from the present (or the past).
It would indicate that if we want evolution equations that involve values just from the past, as in all the heuristic models, we cannot use the framework of the calculus of variations: the motion is no longer along paths that make an action stationary. This is irrelevant for both mathematicians and engineers, but it is irksome for physicists (such a pretty theory. . . ).
But the same can be seen to happen in classical mechanics: as long as the solution is unique for any initial problem at a given time, the solution can be ran backwards or forwards in time. We can express our solution using any time we like: from t = 0 we may predict the values for any t > 0, but from t = t max we can give the values for any t < t max .
But, in that case, we can always use only values from the past to get our solution while in this framework this is no longer possible.
We may change our idea of causality (at least inside the "ideal" world described by mechanics) and allow that, as long as the solution can be determined for all times, it is irrelevant what values are involved: the behaviour of the system is given. This whenever we can establish existence and unicity of the solution which is, in general, still an open problem for many fractional equations.

New Mathematical Scenarios: New Families of Functions and Equations
As mathematical tool, fractional operators establish important relations between transform integrals and special functions. So, the combined use of integral representations, exponential operational rules and special polynomials provides a powerful tool in the formalism of fractional calculus ( [20,21]). Furthermore, fractional operators allow to elude singularities and reduce linear ordinary equations with variable coefficients. As a consequence, an extension of the classical integral representation of the related special functions can be obtained by using fractional operators ( [22,23]).
From a applied point of view, fractional calculus offers a modeling scenario where fundamental mathematical questions converge and appropriate numerical algorithms can de developed. For these reasons, fractional calculus has many applications in different areas [24].
The main contribution of the fractional calculus is the consideration of intermediate-order dimensions through integrals and derivatives of arbitrary order ( [25,26]). This has allowed to get a better modeling in different applications, for instance to model biomedical and biological phenomena ( [27]). A large number of models considering long-range dependence and systems with memory are constructed with integro-differential and fractional equations.
In classical physics, many fundamental equations are based on similar laws: • Hooke's Law: Newton's second law: F(t) = k d 2 x dt 2 (t). As an interpolation of these equations, a fractional approach gives the possibility to look for intermediate or mixed behaviours: Some other contexts are the diffusion processes associated to the basic diffusion equation: as we show in Table 1: Interpolation: Wave equation (hyperbolic): Another fractional approach associated to the previous one is the use of Dirac-type fractional equations [28][29][30].
In this way, the following equation, describes two coupled diffusion processes or a diffusion process with internal degrees of freedom. Depending on the chosen representation of the Pauli algebra, that A and B must verify, we obtain a system of equations coupled or decoupled: In the study of the temporal inversion (t → −t), we have the invariance of the fractional Dirac equation for the values 0 < α < 1: , . . . , 6 7 , 1 9 , . . . Some other fractional differential equations are obtained by considering the root 1/3 of both the wave and diffusion equations: where A possible realization is in terms of the 3 × 3 matrices associated to the Silvester algebra: with ω a cubic root of the unity and Ω a cubic root of the negative unity. In this case, ϕ has three components.

Nonlocal Phenomena in Space and/or Time. Applications
We use the term non-locality if what happens in a spatial point or at a given time depends on an average over an interval that contains that value. Thus, the non-local effects in space correspond to long-range interactions (many spatial scales), while the non-local effects in time suppose memory or delay effects (many temporal scales).
These phenomena are associated to integral or integro-differential equations, which appear in multiple contexts:

•
Potential theory: Newton and Coulomb laws of the inverse of the square of the distance. These different phenomena can be described by fractional differential equations, and it sets out two fundamental questions:

1.
Are the models with space and/or time fractional derivatives consistent with the fundamental laws and symmetries of Nature? 2.
How can the fractional differentiation order be experimentally observed and how does a fractional derivative emerge from models without fractional derivatives?
As example to answer the first question, it is interesting to remark that, for instance, the fractional diffusion equation with some kind of time fractional derivative, verifies the second law of thermodynamics only if the following the generalized Fourier law is satisfied: Not all fractional operators satisfy the condition aforementioned, so it might be the key to choose the convenient fractional derivative or to apply restrictions to the initial and boundary conditions of the problem. Let us define ρ = α − 1. When ρ = 1, the condition Equation (104) is trivial; but when 0 < ρ < 1, this issue is a complex problem with different solutions according to the selected fractional operator and conditions [31].

Application of Fractional Calculus to Model Atmospheric Effects of Absorption
The time-fractional Cauchy problem is well-known: and its solution in the space of functions with Laplace and Fourier transforms, LF = L(R + )xF(R), is defined by where the Mittag-Leffler function is evaluated on the complex plane and G(k) represents the Fourier transform of g(x). For instance, for β = 1 and g(x) = e −µ|x| , µ > 0: The fundamental solution of the problem is obtained for g(x) = δ(x), G(k) = 1 and, in this case, the moments for β = 1 have the following expression: When replacing t by the wave-length of the radiation, λ, the moment for n = 1 and an appropriate constant β returns the Angstrom law, that is used to model the coefficient of molecular scattering τ for the absorption of the incoming energy inside the Martian atmosphere due to the dust [32,33]. The parameters α and β would be fixed in function of the Martian dust features. So, this relation shows a possible application of fractional calculus to model the dynamic of the Martian atmosphere. Deep studies, by using cloud computing, on this issue have been developed in [34][35][36].
In the context of the spatial exploration, other new original application of the fractional calculus analysis is the prediction and identification of dust devils and correlations between wind and seismic signals in a Martian meteorological payload packet. We recently started this project on the basis of our previous experience in the missions to Mars and to apply in ExoMars22 (initially ExoMars20 but now delayed due to the Covid-19) [37][38][39][40].

Chaos in a Fractional Duffing'S Equation
Duffing's equation has been a model for many studies on chaotic systems. It considers a simple but complex system where chaos can appear depending on the values of the parameters. The mathematical model is a time-forced, dissipative, second order nonlinear differential equation that can be viewed as a perturbed Hamiltonian or Lagrangian system [41]. Different possible potentials can be considered but the fundamental equation corresponds to a model for a long and slender vibrating beam set between two permanent magnets, subjected to an external sinusoidal force.
Duffing's equation shows many paradigmatic features of chaotic systems, in a somewhat simple frame, in the Theory of Dynamical Systems. The fractional counterpart we have chosen may possess similar relevant behaviours of other more general fractional models. This is the basis for this study. Besides, although the presence of chaotic behaviour in fractional Duffing's equations has been documented, many questions remain open.
It is possible to extend Duffing's equation into a fractional one in many ways, either as a second order differential equation, or as a system of simultaneous two first order equations, with some or all of these derivatives replaced by fractional ones. Different authors have, thus, considered different fractional equations, playing also with the fractional order of derivation [42,43]. In our case, we have chosen to replace only the first order derivative by a fractional Caputo derivative. The equation we consider isẍ where α ∈ (0, 1) and D α t x stands for the Caputo fractional derivative with lower limit 0 of order α. This equation has the advantage of a regular solution (at least C 2 ) whose existence can be ensured [44]. This is not merely a mathematical model since it can be viewed as the same mechanical device represented by Duffing's equation but immersed in a viscous medium. For some values of the parameters, a strange attractor is obtained, quite similar to the one that appears in the classical (i.e., non fractional) Duffing's equation.
We have studied the controlability of the chaotic regime in the presence of both harmonic and nonharmonic external perturbations, considering geometrical resonances for the second case. Using resonant Jacobi functions, we obtain conditions for external, additional, drivings that ensure chaos-free responses in our model [44].
We have also characterized the chaotic behaviour in our fractional model computing the maximum but, also, all the other Lyapunov characteristic exponents. We have used, as a reference, the fiduciary orbit technique and we have built a perturbative approach with a local equation that allows to estimate all the exponents [45].
The results show that a chaotic regime exists for the fractional Duffing's equation but also a regular regime with a long transient time. This regular regime, in practice, can be assimilated in many cases to a chaotic regime for quite long transient times, although the solutions for even longer times tend to regular, almost-periodic limiting curves.

Conclusions
We deeply thank the enlightening research papers of D. Baleanu and A. Fernandez [4] and M.D. Ortigueira and J.A.T. Machado ( [5][6][7]). These works are a source of inspiration in order to continue the studies about fractional calculus and its applications to real world phenomena. In the light of the results of these works, in this paper we present new scenarios of discussion that might complement the previous studies. The main novelties of this work, can be summarized in the following items: • The effect of the nonlocality, associated to the structure itself of the fractional derivative, manifests that in a solution we can observe the coexistence of two decays exponential and polynomial according to the time scale we consider (Section 1.4).

•
Up to our best knowledge, the obtained results in Section 2 by using the Dirac delta are new. We show that the loss of causality in fractional mechanics is not specific to choosing any particular fractional derivative but is intrinsic to the formulation of nonlocal dependence with more general kernels. We present an approach, using a Dirac delta formulation, that simplifies the, otherwise, more cumbersome computations.

•
Concerning Section 3, there are many remarkable issues. We show explicitly the use of the fractional calculus as an instrument to create new equations as interpolation among other classical ones very well known. An important issue is the interpolation between the parabolic and the hyperbolic dynamics. In this case, we have challenging dynamics attending to the behaviours under discrete symmetries as the time and space inversion.
Dirac obtained his famous equation by considering the square root of the Klein-Gordon equation.
It is related to the basic idea of evolution depending only on the initial configuration of the system. At the same time, Dirac introduced the concept of internal degrees of freedom: the spin of a particle. In this contribution, we apply the above idea of Dirac to the square root of the classical heat equation and we obtain a fractional diffusion equation with internal degrees of freedom.
We extend the idea of considering a general root equation of a given one, and we obtain a connection between the Silvester algebra and the fractional calculus.

•
In Section 4, we show one example where a fractional diffusion equation does not satisfy the second law of Thermodynamics, and we consider the use of the fractional calculus to model the dust dynamics with the associated electromagnetic interaction in Earth and Mars atmospheres.
On the other hand, these developments set out interesting questions and discussions about the adjustment and the reliability of the mathematical models to the dynamic of the real processes; for instance, the objectivity in the descriptions of the models. In this sense, the following items are remarkable: • In physics, the laws must have the same form in all the inertial reference systems, otherwise we could distinguish an inertial reference system from other one by internal experiments.

•
The above statement implies that we must have a suitable dictionary to relate the measurements in one system to other one.

•
For inertial systems we have the Galileo and the Lorentz transforms (dictionaries).

•
Einstein generalized the Galileo transform to the Lorentz transform in order to take into account that two inertial systems cannot distinguish each other by either internal mechanic or electromagnetic experiments (special relativity). • Einstein generalized the above statements to the accelerated reference systems and created the general relativity. • As a consequence, given an evolution fractional differential equation should be analysed its behaviour under Galileo and Lorentz transformations, as well as the discrete symmetries of the time and space inversion. This preliminary analysis will enlighten, for instance, the reversibility and causality issues associated to the equation.

•
Concerning the velocity issue, we consider it in the general sense as the variation of a quantity with time.