Transition from the wave equation to either the heat or the transport equations through fractional differential expressions

We study a model that intermediates among the wave, heat, and transport equations. The approach considers the propagation of initial disturbances in a one-dimensional medium that can vibrate. The medium is nonlinear in such a form that nonlocal differential expressions are required to describe the time-evolution of solutions. Non-locality is modeled with a space-time fractional differential equation of order $1 \leq \alpha \leq 2$ in time, and order $1\leq \beta \leq 2$ in space. We adopt the notion of Caputo for time-derivative and the Riesz pseudo-differential operator for space-derivative. The corresponding Cauchy problem is solved for zero initial velocity and initial disturbance represented by either the Dirac delta or the Gaussian distributions. Well known results for the partial differential equations of wave propagation, diffusion and (modified) transport processes are recovered as particular cases. In addition, regular solutions are found for the partial differential equation that arises from $\alpha=2$ and $\beta=1$. Unlike the above cases, the latter equation permits the presence of nodes in its solutions.


Introduction
Partial differential equations play a relevant role in mathematical physics, mainly those of the second order [1][2][3][4]. In analogy to conics of analytic geometry, they are classified in either elliptic, parabolic or hyperbolic. These equations are useful to describe a diversity of physical phenomena that include wavelike propagation, diffusion and transport processes in practically all the branches of physics. Namely, in continuous and classical mechanics it is common to face problems of either the hyperbolic (vibrating strings, stretched membranes) or the parabolic (heat conduction) types. In statistical mechanics, to describe the motion of a cloud of noninteracting particles, it is necessary to solve a transport (linear Boltzmann) equation, which is of parabolic type, and so on. However, as a fingerprint of the above models, the local character of the differential operators used to construct the corresponding equations permeates over the predictions. The latter produces solutions with unpleasant features from time to time because, occasionally, the predicted behavior is not compatible with the observed phenomenon. For example, the heat equation predicts that thermal disturbances propagate at infinite velocity. Clearly, the model offered by such equation is unrealistic, although it works well when diffusivity is slower than the velocity of propagation. Other examples include the study of systems that exhibit the combination of two (or more) traits. Viscoelasticity, for instance, is understood as the combination of viscous and elastic behavior exhibited by some materials when they are deformed. In such cases the local profile of conventional differential operators represents a limitation to model the involved phenomena in realistic fashion.
A form to remove the intricacies of local operators consists in substituting the usual concept of linear medium (where the phenomenon under study is observed) by the idea of nonlinear medium with memory. This approach has been successfully applied to change the unpleasant infinite propagation velocity of conventional diffusion by a more realistic finite velocity in e.g. [5,6]. The main point is that memory effects may be modeled with nonlocal differential equations that arise from the action of fractional differential operators on the appropriate set of functions. The same approach is useful to investigate viscoelasticity, where the fractional calculus finds a great deal of applications [7].
This work is addressed to study a model that intermediates among the wave, heat, and transport equations. That is, we are interested in connecting pure wavelike propagation with pure diffusion and transport processes by means of the fractional differential properties of a single model in unified form. Our approach considers the propagation of initial disturbances in a one-dimensional medium that can vibrate. When the medium is disturbed at a given time by a vibration, perturbations spread out from the disturbance at a given velocity. The medium is considered nonlinear in such a form that non-locality is required in both, space and time variables. Thus, we are going to study a space-time fractional differential equation with time-derivative of order 1 ≤ α ≤ 2, and space-derivative of order 1 ≤ β ≤ 2. For the time-derivative we shall adopt the one introduced by Caputo [8], the space-derivative will be defined by the Riesz pseudo-differential operator [9]. The values taken by α and β define a squared area in the α × β-plane that will parameterize the solutions of our model. Interestingly, we find regular solutions for the differential equation that arises by fixing α = 2 and β = 1. Unlike the conventional equations aforementioned, this new equation permits the presence of nodes in its solutions. The zeros appear in pairs and propagate, together with the maxima and minima of the solutions, in wavelike form. As far as we know, such equation as well as its solutions have been unexplored in the literature up to now. In this form, the zoology of combined phenomena may be enlarged by involving the above equation and its solutions.
The organization of the paper is as follows. Motivated by the non-locality of fractional differential equations, in Section 2 we formulate a Cauchy problem that embraces, as particular cases, three well known problems of mathematical physics: to solve either the (hyperbolic) wave equation, the (parabolic) heat equation or the (modified) transport equation for properly defined initial conditions. The differential equation associated to α = 2 and β = 1 is a modified version of an equation of the parabolic type which completes the set of second order partial differential equations that can be constructed by omitting the term of simultaneous space and time derivatives. Section 3 includes the solution for two different Cauchy problems, both of them considering zero initial velocity. The first one is defined by the Dirac delta distribution as the initial disturbance. The solutions are written in terms of the H-function, their convergence is analyzed and concrete expressions are given for particular cases that include the well known results of conventional differential equations. Some directrices about the form of intertwining such equations are also given. The second problem is defined by considering a Gaussian distribution as the initial disturbance. We show that the solution is a convergent series of H-functions. To our knowledge, this result has been missing in the literature up to now. We also show that such a series goes to the solution of the Dirac delta case at the appropriate limit. In Section 4 we include the analysis of our results and discuss about the form in which the maxima of the solutions propagate in the t × x-plane. In accordance with other studies, we find that the fractional order of the differential equations affects the behavior of the characteristic integrals in such a form that they are not straight-lines anymore. However, our results show that the time-dependence of the maxima is not as simple as this has been reported by other authors. Some concluding remarks are included in Section 5. To offer fluidity in reading the manuscript, we have resigned detailed calculations of important solutions to appendices B, C, and D. In turn, Appendix A includes useful information and expressions that are recursively used throughout the paper.

Problem formulation
The one-dimensional wave equation [1][2][3][4] with the initial conditions admits a unique solution given by the D'Alembert's formula The Cauchy's problem (1)-(2) is defined for t > 0, with sufficiently smooth functions ϕ 1 and ϕ 2 . The first term on the right-hand side of (3) represents the propagation of the initial wavelike disturbance ϕ 1 , ruled by the characteristic integrals x ± vt = const, for zero initial velocity ϕ 2 = 0. Hence, the positive parameter v represents the velocity of propagation of the wave. Consistently, ϕ 1 (x − vt) propagates to the right and ϕ 1 (x + vt) to the left. In turn, the integral term of (3) corresponds to vibrations produced by the initial velocity ϕ 2 with no initial disturbance ϕ 1 = 0.
A case of special interest is defined by the Dirac delta distribution [3] as follows The initial disturbance (4) is a very localized pulse that levels off rapidly and has a strength equal to µ. Therefore is solution of the one-dimensional wave equation (1) with the initial conditions (4).
On the other hand, to define uniquely the solution of the one-dimensional heat equation it is generally sufficient the initial condition [10] u(x, 0) = ϕ 1 (x), where u(x, t) is the absolute temperature and k > 0 stands for the thermometric conductivity (or diffusivity) [10]. Any line parallel to the x-axis is a characteristic integral of Eq. (6). The Cauchy problem (6)- (7) is useful to describe diffusion processes without drift in one-dimension. If the initial function (7) is the Dirac delta distribution then (6) is the Fokker-Planck (or forward Kolmogorov) equation [11]. Indeed, it may be shown that the time-dependent Gaussian density often called source solution [1], satisfies (6) with the initial condition Clearly, density (8) converges to the initial distribution (9) at the limit t → 0 + . Explicitly Besides, at fixed time, u S (x, t) falls off very rapidly as |x| increases. Therefore, density (8) is useful to describe diffusion in one-dimensional media where the heat is concentrated in the vicinity of x = 0 [4,10]. In addition, u S (x, t) is positive and nonzero for any x and t, so one also finds that "the effect of introducing a quantity of heat at x = 0 is instantaneously noticeable at remote points" [10]. That is, thermal disturbances propagate at infinite velocity. Although unrealistic, this model works well in situations where diffusivity is slower than the speed of propagation. If diffusivity and the speed of thermal waves are compatible then (6) must be modified to include an additional second order time-derivative term [3].
One of the directions of this work is addressed to study a model that intermediates between the wave and the heat equations described above. That is, we look for a connection between pure wave propagation and pure diffusive phenomena in one-dimension. A first resource consists of replacing the second order time-derivative of the wave equation (1) by a fractional time-derivative D α of order α ∈ [1,2]. In this case we have Here, v α > 0 is a constant expressed in the units standing for length and time units. Using α = 2 in (11) we recover the wave equation (1) while α = 1 gives the heat equation (6), with v α=2 = v and v 2 α=1 = k respectively. Remarkably, the unpleasant feature of the heat equation (6) that thermal disturbances propagate at infinite velocity is removed since (11) describes heat conduction for nonlinear materials with memory, where the speed of propagation is finite [5,6]. Extending the values of α to include the sub-diffusion interval 0 < α < 1, one arrives at a fractional version of the Fokker-Planck equation for the initial condition (9). However, such a case is out of the scope of the present work and will be studied elsewhere. Preliminary results associated to (11) have been already reported by other authors in e.g. [12][13][14][15][16][17] (see also the general approaches discussed in [7,11,[18][19][20][21]).
To include a major diversity of cases we may replace the second order space-derivative of the heat equation (6) by a fractional space-derivative D β of order β ∈ [1, 2] to obtain The constant k β > 0, measured in units [k β ] = [L] β [T ] −1 , is such that k 2 = k. Equation (12) is useful to describe symmetric Lévy-Feller diffusion processes in which jump components are present [11,22] (see also [23,24], and general approaches in [7,11,[18][19][20][21]). The parameter β = 2 gives the heat equation and β = 1 leads to a modified version of the one-dimensional transport equation [25,26] It may be shown that for any smooth function f (x), the expression u(x, t) = f (x + k 1 t) solves the conventional transport equation; therefore, for the initial condition (7) we have u(x, t) = ϕ 1 (x + k 1 t), where k 1 is the fluid (hydrodynamic) velocity [27]. The latter is not automatic for Eq. (13), see [25,26]. For the sake of simplicity, Eq. (13) will be referred to as transport equation.
The most general alternative corresponds to substitute both derivatives, the temporal and the spatial ones, by their fractional versions in Eq. (1). Thus, we have the space-time fractional differential equation [25,26] where the coupling constant v α,β > 0 is written in the units [v α,β ] = [L] β/2 [T ] −α/2 , and is such that v α,2 ≡ v α , v 2 1,β = k β , and v 2 2,1 = κ. The above equation can be studied in different approaches [7,11,[18][19][20][21], and offers very interesting mathematical features as the fact that one can transit from the wave equation (1) to the heat one (6), and vice versa, by fixing β = 2 and running α in [1,2], see segment BA in Figure 1. In turn, the heat equation (6) serves as point of departure to arrive at the transport equation (13) by fixing α = 1 and permitting β to take values in [1,2], see segment AD in Figure 1. The transition from the transport equation (13) to the wave equation (1) is reached by making α = β, with α ∈ [1, 2], see segment DB in Figure 1. Of course, any other (differentiable) path described by the point (α, β) to connect B with either A or D is admissible in the α × β-plane. As indicated above, we are interested in those paths described within the grey squared area of Figure 1, and the segment-lines I-V are representative of the simplest ones. In this respect, the zone enclosed by the circuit BAD refers to phenomena that intermediate among wave propagation, diffusive and transport processes. The zone delimited by the circuit BCD represents an additional resource of information. We shall discuss the subject in the next sections. Figure 1: The α × β-plane associated to the space-time differential equation (14). The square 1 ≤ α, β ≤ 2 contains the points (α, β) that are the subject of interest in this work. Vertices B, A and D correspond to the wave equation (1), the heat equation (6), and the transport equation (13), respectively. The solutions of (14) for different initial conditions (2) associated to the grey squared area are defined within the text. Special attention is paid to the five segment-lines defined by α = 1, 2, β = 1, 2, α = β, and the vertices A, B, C, D. See complementary information in Table 1.
In addition to the differential equations (1) and (6), together with Eq. (13), we remark the differential form which is defined by the space-time fractional differential equation (14) evaluated at vertex C, and will be called complementary equation. This can be achieved from either the wave or the transporttransport equations, (1) and (13), through segment-lines II and III, respectively. In this form, the circuit ABCD is completed as the combination of triangles BAD and BCD of Figure 1.
As we can see, the diversity of points (α, β) that can be used to define a concrete form of Eq. (14) embraces important cases of the family of second-order partial differential equations where the sub-labels of u refer to partial derivatives in conventional notation. The term "complementary" coined for Eq. (15) is then justified. In the language of linear second order partial differential equations the set (16) includes only hyperbolic and parabolic types. Indeed, the wave equation (1) is hyperbolic while the heat equation (6) is parabolic. We shall say that complementary and transport equations, (15) and (13) respectively, are also "parabolic".
For the time-derivative D α we shall adopt the one introduced by Caputo [8]. In turn, the Riesz pseudo-differential operator D β [9] is going to substitute the space-derivatives in our approach.

Solution and examples
In this section we study the Cauchy problem defined by the space-time fractional differential equation (14) and the initial conditions (2). In general, we assume that the solutions u(x, t; α, β) : R × (0, ∞) → R will represent the propagation of initial wavelike disturbances ϕ 1 for the initial velocity ϕ 2 in a medium that can vibrate. That is, when the medium is disturbed at a given (initial) time by a vibration ϕ 1 , perturbations spread out from the disturbance at initial velocity ϕ 2 . From now on we omit the explicit dependence of u on the parameters α and β for the sake of simplicity. It is expected that u will converge either to the solutions of the wave equation (1) or to the solutions of the heat equation (6) at the limits (α, β) → (2, 2) and (α, β) → (1, 2), respectively. The same holds at the limit (α, β) → (1, 1) for the transport equation (13). Our main interest is addressed to the values 1 < α, β < 2, with either α = β or α = β. In Figure 1 we show the α × β-plane associated to the functions u(x, t) we are interested in.
Let us start by calculating the Laplace transform L of Eq. (14) in the time-variable. After introducing the initial conditions (2) we arrive at the new equation where U (x, s) stands for the Laplace transform of u(x, t). Now, we calculate the Fourier transform F of Eq. (17) in the space-variable. We obtain which is immediately solved by the function In the above equations we have introduced the notation To retrieve U we calculate the inverse Fourier transform Providing the initial disturbance ϕ 1 and velocity ϕ 2 , the solution we are looking for is obtained from the inverse Laplace transform Next, we determine the explicit form of u(x, t) for ϕ 2 = 0 and two different functions ϕ 1 .

Delta-like disturbances
For initial conditions (4), the integral (20) is simplified as follows The straightforward calculation (see details in Appendix B) gives where H m,n p,q [z|−] ≡ H(z) stands for the Fox function, H-function for short [28,29] (see also Appendix A). The absolute value of x in (22) obeys parity properties of the Fourier transform. The full derivation of the inverse Laplace transform of (22) is included in Appendix C, we obtain Equation (23) represents the propagation of the disturbance produced by a very localized pulse ϕ (δ) 1 (x) in one-dimensional media, according to the space-time fractional differential equation (14).

Discussion about the convergence of solutions
To analyze the convergence of the functions u(x, t) defined in (23) one may consider the series expansion of the H-function given in Theorems 1.3 and 1.4 of Ref. [28] (see Appendix A). It can be distinguished two cases. For α > β (points delimited by the triangle BCD of Figure 1) we have the series which converges absolutely [25]. On the other hand, for β > α (points bordered by the triangle BAD of Figure 1), the function u(x, t) acquires the form Thus, any point (α, β) in the zone defined by the circuit BCD produces a convergent solution (23). In turn, for the points (α, β) within the triangle BAD, the convergence of u(x, t) must be determined from the series (25). For instance, as α = 1 produces the cancelation of the first term on the right-hand side of (25), we arrive at the following expression for the points on segment-line IV of Figure 1, Remark that, in agreement with the results reported in [25], Eq. (26) do not include the point (α, β) = (1, 1) associated to vertex D.
The case α = β (segment-line V of Figure 1) intermediates between the two series aforementioned. The straightforward calculation [25] shows that u(x, t) is reduced to the expression

Special cases
We find the following special cases of function (23). The sub-label of u refers to the diagram shown in Figure 1. For solutions along the segment-lines I-V we have where φ(a, b; z) is the Wright function (see Appendix A). The above expression is a well known result for Eq. (11) with initial conditions (4), see e.g. [25]. We also have According to Section 3.1.1, the above solutions are absolutely convergent series. However, although u II and u III give rise to u B and u D respectively, they are scarcely explored in the literature. On the other hand, both of them converge to the solution associated to vertex C. That is, they lead to the solution of the complementary equation (15) which, in turn, seems to be missing in the set of exactly solvable differential equations that are commonly studied in physics and enginery (see for instance the monograph [21]).
The function converges to u A (x, t) and u D (x, t) as β → 2 and β → 1, respectively. As solution of the space-differential equation (12), the expression (31) is in correspondence with results that have been already reported in e.g. [11,[22][23][24]. The comparison of (31) with (26) gives rise to the following evaluation of the H-function The function embraces the results for α = β reported in e.g. [25,26], and regulates the simplest transition from u D (x, t) to u B (x, t) and vice versa. Comparing (32) with (27) we obtain another evaluation of the H-function. Namely |z| β−1 γ β sin(πβ/2) γ 2β + 2|z| β γ β cos(πβ/2) + |z| β , γ > 0. (33) The solutions defined by the four vertices can be written as Compare Eq. (34) with u I (x, t) and u IV (x, t). The function can be recovered from either u I (x, t), u II (x, t) or u V (x, t). For the functions and see u II , u III , and u III , u IV and u V , respectively. In Table 1 we summarize the above results.

Recovering conventional results
Wave equation. To evaluate u(x, t) at the limit (α, β) → (2, 2) we use the solution (35) associated to vertex B. Of course, u I , u II , and u V are also useful at the appropriate limits. Using Eqs. (A-4) and (A-5) of Appendix A we get

From (A-3) of Appendix A we see that the above integral corresponds to the inverse Mellin transform of µδ(|x| − vt)/2. Then
This last expression reproduces the solution u δ (x, t) given in (5) for the one dimensional wave equation (1) with initial conditions (4).
• As collateral result, the combination of (35) and (39) gives rise to the following evaluation of the H-function: Heat equation. In similar form, the limit (α, β) → (1, 2) is obtained from the solution (34) associated to vertex A. Equivalent results are obtained from either u I or u IV at the appropriate limits. We have Using the duplication formula of the Gamma function (A-17), as well as the Cahen-Mellin integral (A-18), we arrive at the function u S (x, t) defined in (8), which is solution of the one-dimensional heat equation (6) with initial condition (9). The verification that u A (x, t) satisfies the initial condition (9) is easily achieved by using Eq. (10).
• In this case, the combination of (34) and (8) gives rise to the evaluation From (10) we also have Transport equation. The limit (α, β) → (1, 1) is evaluated from (37), associated to vertex D (other options include either u III or u IV at the appropriate limit). That is, from (A-4) and (A-5) of Appendix A, after using the duplication formula (A-17), we have From the Mellin transformation (B-3) of Appendix B, one gets Then, (44) is reduced to the curve which is consistent with Eq. (27). The latter result is well known in the literature (see, e.g. [25], Remark 3.2). In the space variable, u D (x, t) describes a bell-shaped curve known as either Cauchy (mathematics), Lorentz (statistical physics) or Fock-Breit-Wigner (nuclear and particle physics) distribution [30]. It is centered at x = 0 (the location parameter), with a half-width at half-maximum equal to k 1 t (the scale parameter) and amplitude (height) equal to µ(πk 1 t) −1 . That is, at time t, the fluid velocity k 1 defines the width of the disturbance between the half-maximum points x = ±k 1 t. The function u D (x, t) satisfies the Dirac delta constraints (10), which can be verified at the elementary level. Therefore, (45) satisfies the transport equation (13) with initial condition (9).
• The combination of (37) with (45) gives the evaluation with the following limit Complementary equation. For α = 2 and β = 1 the series (24) is simplified as follows Depending on the value of ξ = t 2 /|x| > 0, the alternation of sign in Γ(−n− 1 2 ), n = 0, 1, . . ., produces pairs of zeros in u C (x, t) that are located symmetrically with respect to x = 0. As an immediate example, for ξ << 1 the series (48) coincides with the peaked density which is consistent with the initial condition ϕ 1 (x). The number of zeros increases as ξ → ∞. We will find a first pair at second order of ξ 2 , and so on. On the other hand, u C (x, t) goes to zero from above as |x| → ∞, see Eqs. (48)-(49). So that the zeros produced during the propagation of the pulse ϕ

Other representations
It is well known that the H-function embraces a large number of functions of the hypergeometric type [29]. Among the functions preferred in fractional calculus, besides the H-function, one finds Wright and Mittag-Leffler functions [7,[18][19][20]. To show the flexibility of our results, in this section we rewrite the function (23) in terms of the generalized Wright function. The translation to other hypergeometric representations is straightforward.
After changing z by 1 − s, and using the duplication formula of the Gamma function (A-17) in the Mellin-Barnes integral representation (C-1) of u(x, t), we have The change of variable s = βy produces (51)

Using Eqs. (A-5)-(A-4) of Appendix A one can write
Now, from Eq. (A-16) of Appendix A we finally arrive at the function . (53) The Mellin-Barnes representation of u(x, t) given in (50) coincides with the expression of the Green function reported in [25] for the Cauchy problem defined by Eq. (14) and the initial conditions u(x, 0) = ϕ 1 , ϕ 2 = 0. In agreement with our results, in [25] it is found that the concrete form of u(x, t) given in (53) is reduced to the absolutely convergent series (24) for α > β. Other values of α and β lead to the expressions discussed in Section 3.1.1 and Ref. [25].

Gaussian disturbances
The Gaussian density (8) offers a very versatile profile to define the initial condition ϕ 1 . Namely, at arbitrary time t = t 0 > 0, we can take What makes density (54) interesting as initial condition is that, unlike the Dirac delta distribution ϕ (δ) 1 (x), this is finite at x = 0 for x 0 = 0. The configuration studied in the previous sections for ϕ (δ) 1 (x) can be recovered from (54) at the limit x 0 → 0 + , as it is clear from Eq. (10). Additionally, we have shown that the Gaussian density (54), as well as the Dirac delta pulse (9), can be expressed in terms of the H-function, Eqs. (42) and (43) respectively. The full derivation of the solution of Eq. (14) with initial conditions (54) and ϕ 2 = 0 is given in Appendix D. The final result can be written as the series where Remarkably, with exception of the term with k = 0, the coefficients of the series (55)-(56) become zero at the limit x 0 → 0 + . Therefore we arrive at the expression which coincides with the solution (23) of the space-time fractional differential equation (14) with initial conditions (4), as expected.

Analysis of the results
As we have indicated in the previous section, using the Gaussian density (54) as initial condition is, in many respects, more realistic than using the Dirac delta distribution (9). The main point is that ϕ 1 (x) diverges as |x| → 0. However, the disadvantage of ϕ 1 (x) is that the former, although centered at x = 0, spreads along the entire x-axis while the latter is very localized at x = 0 and equal to zero for any x = 0. Nevertheless, ϕ (δ) 1 (x) may be recovered as a Gaussian density with infinitesimally-narrow width. In the panel of Figure 2 we show the propagation of the initial disturbance ϕ  is closest to vertex A. The clearest feature of the wave-propagation is the decoupling of the initial disturbance into two perturbations, the maxima of which evolve in time according to the respective characteristic integrals. As indicated above, these integrals are the constants x ± vt for vertex B. In our case, the characteristic integrals are not straight-lines (see below for details). In turn, diffusion is characterized by the fading of the initial disturbance as time goes pass (conservation principles imply that the disturbance must spread out with time in order to preserve the area under the initial curve ϕ 1 ). The behavior of the functions depicted in figures (b) and (e) shows, with clarity, a mixture of diffusion and wave propagation.
A second group, formed by figures (c-i) and (e), shows the split of the initial disturbance into two perturbations that is characteristic of the wave-propagation. Additionally, the perturbations take negative values producing the presence of nodes as they propagate. The latter is markedly notorious in figure (i), where (α, β) is closest to vertex C. As we Figure 3: Propagation of the (right hand side) maxima of the solutions to the time-fractional differential equation (11) with initial conditions (54) and ϕ 2 = 0. We have taken u e (x, t) defined in (55)-(56) with β = 2, see Figure 2. The plots correspond to α = 1.9 (solid square, black), α = 1.7 (disk, blue) and α = 1.5 (solid triangle, red).
have shown, such behavior is a trait of the solutions of the complementary equation (15). Therefore, the behavior of functions included in this group is a mixture of wave-and complementary-propagation.
We have two additional groups, the behavior of which exhibits respectively a mixture of transport and complementary-propagation, and a mixture of diffusion and transport processes. The former group includes figures (i-g) and (e), and the latter is integrated by figures (g-a) and (e). With respect to the characteristic integrals, in Figure 3 we show the distribution of the maxima x M (t) of u e (x, t) in the plane t × x for β = 2 and three different values of α. That is, Figure 3 refers to the propagation of the Gaussian density ϕ (G) 1 (x) defined by the points (α, β) along the segment-line I. As we can see, for α = 1.9 the characteristic integral defines a path that is almost a straight-line, which might be expected since (α, β) = (1.9, 2) is very close to vertex B, which corresponds to the wave equation (see the pretty explanation of the behavior of characteristic integrals for wave equation in [4]). On the other hand, the path defined by (α, β) = (1.5, 2) is far away from a straightline. The results depicted in Figure 3 have been calculated numerically, and obey a rule (obtained by the best-fit technique) of the form where c α and θ α are parameterized by α. In contrast with other works like [14], we have found that θ α = 0 and c α = 1 for 1 < α < 2. Of course, c 1 = 0 and c 2 = 1, with θ 2 = 0.
On the other hand, in Figure 4 we show the behavior of the maxima x M (t) for the solutions defined along the segment-line V . In this case, we have found the rule with c α and τ α determined by α. The value of x M (t) is very close to that of x M (t) for α = 1.9, as expected. As in the previous case, unlike other works [14], in our approach the rule obeyed by the maxima is not as simple as t α/2 . To conclude our analysis, the panel of Figure 5 shows the behavior of (23) for six different points (α, β). They correspond, in clockwise direction, to the functions depicted in Figures 2(a)-2(c), 2(e), 2(g), and 2(d), respectively. Here, Figure 5(f) is the barycenter. Unlike the Gaussian case, in the group exhibiting propagation-diffusion behavior, Figures  5(a)-5(c) and 5(f), the wavelike propagation is markedly differentiated at short times. Notice that the transport-diffusion processes are also "accelerated" in the sense that changes are presented in intervals of time that are shorter than those spent by the gaussian profiles. A similar phenomenon is observed for the other groups (not included in Figure 5).

Concluding remarks
The discussion of Section 4 motivates the addition of a new segment-line in Figure 1, named V I, connecting the vertices A and C. Thus, V I defines the simplest path to transit from the heat equation (6) to the complementary equation (15). This can be written as the rule β = 3 − α for 1 ≤ α ≤ 2. The presence of the new vertex E splits the square into four different triangular areas, respectively bordered by the circuits ABE, BCE, CDE, and DAE, see Figure 6. These regions of the α × β plane are in correspondence with the groups mentioned in Section 4, so they serve to classify the solutions of the space-time fractional differential equation (14). Namely, solutions associated to points within the region ABE will exhibit combined properties of wavelike propagation and diffusion. Those associated to BCE will show a mixture of wavelike and complementary propagation, and so on. We have shown that the H-function permits to express, in unified form, a family of solutions to the Cauchy problem defined by the space-time fractional differential equation (14) and the initial conditions (2). We have specialized to the case ϕ 2 = 0 for two different forms of ϕ 1 . To be concrete, we have solved the problem for the Dirac delta pulse (9) and for the Gaussian density (54) as well. We have shown that the solutions of the former can be obtained from those of the latter at the appropriate limit. In contrast with the Dirac delta distribution, which is recurrently used as initial condition in different fractional approaches, the Gaussian density is rarely included as initial condition in the literature on the matter. The material included in this work is addressed to fill such lack of information.
The well known solutions of three important equations studied in mathematical physics have been recovered as particular cases. We refer to the wave equation (1), the heat equation (6), and the transport equation (13), which respectively correspond to vertices B, A, and D of Figure 1, see also Figure 6. In addition, we have completed the set of conic-type second-order partial differential equations F (x, t, u, u x , u t , u xx , u tt ) = 0 by including the differential equation that results from making α = 2 and β = 1 in the space-time fractional differential equation (14). That is, we have solved also the differential equation associated to vertex C, the explicit form of which is introduced in (15) and has been called complementary equation throughout the present work. In this form, besides including the fractional cases where (α, β) does not coincide with any vertex, our approach unifies in a single expression the solutions of four different problems that are usually treated in separate form in mathematical-physics.
Collateral results include the evaluation of different forms of the H-function. As far as we know, most of them have been unclassified up to now in the specialized literature. Particularly, the expressions of the H-functions in which the Dirac delta distribution is involved. We hope these results will be useful for researchers in the area.
The approach presented here may be extended to include nontrivial initial velocities ϕ 2 = 0. Specifically, the Dirac delta distribution ϕ (δ) 2 (x) = γδ(x − x 0 ) offers the possibility to calculate the Green functions for the problems above mentioned. Namely, it is well known that the calculation of the Green function for a given Cauchy problem is equivalent to the calculation of the solutions to the homogeneous equation with the boundary conditions ϕ 1 (x) = 0 and ϕ (δ) 2 (x) [3]. The latter works well for vertices A, B and D, which produce conventional partial differential equations. The verification of such a property for vertex C and any other point (α, β) in the squared area 1 ≤ α, β ≤ 2 is an open problem that we shall face elsewhere.
To conclude this work we would like to emphasize that the substitution of conventional derivatives by their fractional versions in a given dynamical law produces the emergence of interactions that are not apparent (and cannot be noticed) in conventional models. The framework developed here look for connections between different laws that are already known. But fractional calculus offers a more wide range of possibilities. For example, in the case of the harmonic oscillator, the fractional time-derivative adds a sort of frictional forces that are not justified in the conventional Newtonian model since no environmental interactions are assumed a priori [31]. Therefore, the fractional version of the oscillator presupposes either that the system suffers a kind of self-interaction or that it is embedded into a medium with memory. Both assumptions predict phenomena that may require a new dynamical law for their explanation. A similar situation occurs for the fractional quantum oscillator [32], for which there exist some immediate applications in quantum optics [33]. Work in this direction is in progress.

A Useful definitions and expressions
• The Mellin transform [34] of a function f (x) is defined as

Consistently, the Mellin-Barnes integral (inverse Mellin transformation) is written as follows
In particular, for  A 1 ), . . . , (a p , A p ) where m, n, p, and q are nonnegative integers, a i , b j ∈ C, A i , B j ∈ (0, ∞), and The contour L in (A-4) separates the poles of Γ(b j − B j s), j = 1, . . . , m, from the poles of Γ(1 − a i + A i s), i = 1, . . . , n. Detailed information can be found in [29].

C Derivation of Eq. (23)
The inverse Laplace transform of U (x, s) is calculated at the elementary level by noticing that the s-variable, which is to be transformed, is encoded in the term λ 1−z /s appearing as a factor on the right hand side of (B-6). From (B-2) we have