Root Operators and “Evolution” Equations

Root-operator factorization a la Dirac provides an effective tool to deal with equations, which are not of evolution type, or are ruled by fractional differential operators, thus eventually yielding evolution-like equations although for a multicomponent vector. We will review the method along with its extension to root operators of degree higher than two. Also, we will show the results obtained by the Dirac-method as well as results from other methods, specifically in connection with evolution-like equations ruled by square-root operators, that we will address to as relativistic evolution equations.


Introduction
Evolution equation is a rather generic term, which in Physics signifies any mathematical device useful to describe the evolution of a dynamical system.Here, we will adopt this term to refer to a partial differential equation of the type [1] for the dependent variable V, designated to characterize the physical state of the dynamical system of concern.Equation (1) specifies in fact the rate of change of V, regarded as a function of the independent "space" and "time" variables (x, t), with respect to the relevant evolution variable "t".
The left hand side is just the first order derivative of V with respect to t, whilst the right hand side, which may be both linear and non linear, only involves V, its integer-order derivatives with respect to x, and possibly the independent varables x, t.
The meaning of the "state function" V(x, t) depends on the specific problem at hand as also the definite form of the right hand side.Anyway, since the equation only involves the first order time derivative, one expects its solution to be uniquely specified by a single initial condition V(x, t 0 ) = V 0 (x), so that, known V at a given time t 0 , it can be known at any subsequent time t > t 0 after integrating the equation.
Many physical processes are well modelled by equations like Equation (1).The heat equation (HE) [2], the time-dependent Schrödinger equation (SE) [3] and the paraxial wave equation (PWE) [4] are substantive examples of evolution equation.The theory concerned with such equations is well established; we will review some basic aspects of it below, thus exemplifying the essential features of the theory concerned in general with evolution-type equations.
Many physical processes are equally well modelled by equations, which cannot be traced back to the scheme of Equation (1).In fact, the HE, SE and PWE pertain to specific physical contexts and also to definite approximations within those contexts.Thus, the HE pertains to the non-relativistic theory of heat diffusion; the SE pertains to non-relativistic quantum mechanics as the PWE to the paraxial wave optics.
When trying to go beyond the inherent approximations in order to improve the adherence of the equation to the process it describes, we may be led to face equations involving higher-order derivatives of V with respect to the "time" variable t, or equations involving fractional differential operators (modern theories of transport in heterogeneous porous media resort, for instance, to fractional advection-diffusion equations) or pseudo-differential operators.
The relativistic heat equation (RHE), as proposed in the telegraph form in [5], α and C being the thermal diffusivity and the spead of heat, the Klein-Gordon equation (KGE) [6] involves the square-root of an operator containing the laplacian ∇ 2 .As is well known, Equation ( 5) is also referred to as Salpeter equation [8]; it has been the object of recent investigations (see, for instance, [9][10][11]).
The analysis here presented is only an aspect of an investigation, which resorting to different methods, is aimed at establishing if, at which extent and in which form some properties of the evolution equations could be recovered to equations, which are not of evolution type or demand to deal with pseudo-differential operators, like the aforementioned RHE, KGE,WE and RSE.
As shown in [12][13][14], the Dirac-like factorization approach conveys a valuable method to tackle with both kinds of difficulties.
In fact, the Dirac equation [6,15] γ j ∂ ∂x j − mc h I 4×4 ψ = 0, j = 0, 1, 2, 3 where γ j , j = 0, 1, 2, 3, are the 4 × 4 Dirac matrices, (x 0 , x 1 , x 2 , x 3 ) ≡ (ct, x, y, z) and I 4×4 is the 4 × 4 unit matrix, offers in a sense the "evolution-like" alternative to the KGE.As is well known, it was originally formulated by Dirac when seeking a relativistically covariant evolution equation for the state function of a quantum particle in the Schrödinger-like form, i.e., of the type with the Hamiltonian H being a linear Hermitian operator [15].However, it can also be understood as following from a "factorization" of the Klein-Gordon equation so that one can eventually deal with a first-order derivative with respect to the evolution variable even though in a system of four coupled linear differential equations for the 4-component state vector ψ.The Dirac-like factorization approach can as well be effectively applied to deal with evolution equations ruled by fractional differential operators or pseudo-differential operators, like that entering the RSE.In fact, allowing to express the power of operators as the sum of operators, it allows for the "disentanglement" of root operators into the sum of operators, and hence, under appropriate conditions, one can overcome the problem of working with fractional differential operators [12].
In addition, the factorization approach to root operators may open new perspectives within the theory of fractional calculus, suggesting, for instance, alternative formulations to already well-established definitions and/or treatments [13].
The plan of the paper is as follows.In Section 2 we will briefly recall the main features of the theory concerned with evolution equations, referring to the HE, the SE and the PWE as basic examples.Section 3 synthesises the relation between the evolution operator formalism and the solution of fractional partial differential equations.In Section 4, we will review the Dirac-like factorization method in connection with root operators.Square, cube and quartic root operators will be treated in some detail, specific properties of the relevant Lie algebra of inherent matrices will be deduced.Then, we will illustrate some applications of the method, specifically in connection with the square and cube root operators in Section 5. Finally, in Section 6 we will focus the discussion on evolution equations ruled by relativistic-type free Hamiltonian operators, that we will address to as relativistic-type free evolution equations.We will show that the solution to such equations can be expressed in the form of an integral transform of the initial data in full analogy with the HE, which can in a sense be considered as the "non-relativistic" counterpart.Some properties of the obtained solution will be deduced tracing a comparison with those of the solution of the HE.The concluding notes of Section 7 will close the paper.

Examples of Evolution Equation: Heat, Schrödinger and Paraxial Wave Equations
The (1+1)D HE [2] rules the time evolution of the temperature function u(x, t) through an in principle infinitely long homogeneous bar, characterized by the thermal diffusivity α.The relevance of such an equation is not limited to the context of heat propagation, and not only to the context of second order linear evolution equations.For instance, the non linear Burgers' equation [16] can be converted into the HE by the Hopf-Cole transformation [17,18].Also, the Black-Scholes equation [19,20] can be trasformed into the HE; originally proposed in the early 1970's as a model for investment portfolios, the Black-Scholes equation is now at the heart of the modern financial industry.At a basic level, we write down the (1+1)D SE [3] which, as evolution equation for the wavefunction ψ(x, t), describes the 1D free motion of a (spinless) particle of mass m.Correspondingly, the 2D PWE [4,21] is the equation of motion for the complex slowly-varying amplitude ψ(x, z) of a monochromatic scalar light field propagating in a homogeneous medium of refractive index n.Here, z denotes the coordinate along the main direction of the field propagation, k 0 = 2π/λ 0 is the field wave number in the vacuum, and ω 0 the relevant angular frequency, ω 0 = k 0 c.The analogies between the two equations are evident, and may be synthesized by the correspondences Likewise, the analogy with the HE becomes recognizable once allowing the evolution variable in Equations ( 8) and ( 9), or equivalently in Equation (7), to take on a purely imaginary value.
The theory concerned with the heat, Schrödinger and paraxial wave equations is well established.However, in spite of the afore-evidenced analogies, the three equations have been analysed by different approaches, yielding parallel formulations, that only recently have been merged into each other within the framework of a united formalism.

General Treatment: Hamiltonian and Evolution Operators
Evidently, the HE, the SE and the PWE can be recast in the form introducing a sort of "Hamiltonian" operator H (not necessarily Hermitian), which can be seen as the "generator" of the system evolution.Specifically, the "Hamiltonian" pertaining to the aforementioned equations is simply the free-evolution operator, The solution to Equation ( 10) can be written in the form introducing the "evolution operator" U (not necessarily unitary), which acts on the initial data ϕ(x, 0) = ϕ 0 (x) to yield the wavefunction at subsequent times.Specifically, for the equations we are considering the evolution operator is and the "time-like" variable ζ is intended to be ζ = αt, ζ = iht 2m and ζ = iz 2k 0 n according to the equation under consideration.

Evolution Operator as Poisson and Fresnel Transforms
Under the minimal assumption that the initial condition ϕ 0 (x) tends to zero sufficiently rapidly as x → ±∞, the evolution operator specializes into the Poisson transform for the HE and into the Fresnel transform for the SE and PWE [22].
Thus, the solutions to the (1+1)D HE follow as [2]: where u 0 (x) = u(x, 0) conveys the initial condition.The kernel of the transform is the fundamental solution of the HE, i.e., the solution corresponding to a "point-like source" at t = 0: u 0 (x) = δ(x), signifying a highly concentrated unit heat source; in fact, lim t→0 S(x, t) = δ(x).Analogously, the solutions to the (1+1)D SE and the 2D PWE result from the Fresnel transform of the initial condition ψ 0 (x) = ψ(x, 0) as with ς being ς = ht 2m or ς = z 2k 0 n according to whether we are interested in Equations ( 8) or (9).The above can be seen as a sort of "imaginary" version of Equation (15).Again, the kernel is the fundamental solution of the equations of concern, being, as before, It is worth noting that Equations ( 15) and ( 17) implicitly assume the equivalence of the representations of the evolution operator U as an exponential operator (Equation ( 14)), involving the free-Hamiltonian operator ∂ 2 ∂x 2 , and as a Poisson or Fresnel transform.Evidently, the former requires the implied series of derivatives of the initial data, ∂ 2n u 0 ∂x 2n or ∂ 2n ψ 0 ∂x 2n , to exist and converge to a finite value, whereas the latter is meaningful only if u 0 (x) or ψ 0 (x) are integrable and the pertinent integrals converge.The discussion of the legitimacy of such an equivalence is out of the topic of the paper [22].

Polynomial Solutions and Symmetry Transformations
In view of the considerations we will develop in Section 4 in connection with the treatment of relativistic-like evolution equations, we will focus here on two aspects of the theory of the HE, the SE and the PWE: The polynomial solutions and the symmetry transformations.
The polynomial solutions v n (x, t) of the HE, referred to as heat polynomials [2], correspond to initial data given by monomials: v n (x, 0) = x n , and hence For simplicity's sake, we have set α[m 2 /s] = 1.Also, an appropriate multiplying constant of unit value is implied in the initial monomials in order to provide the v n s with the correct dimensions in conformity to those of u(x, t).
The heat polynomials realize the power series expansion of the simple exponential solution of the (1+1)D HE [2] with χ arbitrary parameter.In practice, one considers the eigenstate e χx of the derivative operator ∂ x belonging to the eigenvalue χ.It evolves indeed according to Equation (20), also obtainable from the inherent power series expansion, by evolving the sinlge monomials, i.e., Indeed, a significant aspect of the theory of the HE is aimed at establishing the criteria under which the power series expansion of a function (when it exists) can be used to determine the evolution of that function, when given as input for the HE.
Equation (20) yields for the v n 's the explicit expression with H n denoting Hermite polynomials [26].
For future use, let us write down the explicit expressions of the first v n s: A further property of the v n s of interest here concerns the raising and lowering operators (or, multiplication and derivative operators), conveyed by the relations yielding the differential equation obeyed by the v n s, Recently, polynomial solutions of the 2D PWE have been introduced in full analogy with the heat polynomials [27,28].In fact, in the light of the afore-mentioned analogy, the optical analogues v n (x, ς) of the heat polynomials are found to be explicitly given by Also, they result from the Fresnel transform of monomials, i.e., y n dy, and realize the power series expansion of the function which simply conveys the plane-wave solution of both Equations ( 8) and ( 9), corresponding to the input E(x, 0, λ) = e iλx , i.e., the eigenstate of the momentum operator −i∂ x , belonging to λ.The parameter λ is then related to the transverse wavenumber (i.e., spatial frequency) of the field or to the momentum of the particle; the frequency chirping in optics corresponds indeed to the "energy chirping" in quantum mechanics.The investigation of the symmetry transformations, pertaining to the equations we are dealing with, is a very fascinating and fruitful job.Symmetry transformations establish specific rules to pass from one solution to another solution of the same equation.
In particular, in relation with the HE, SE and PWE, we may quote the special conformal transformation that allows for an appropriate Gaussian modulation of a solution to get another solution.In fact, it can easily be verified that if u(x, t) is a temperature function, also is a temperature function for any arbitrary parameter ε.It can be traced back to the Gaussian modulation of the initial data, which then turn from u 0 (x) to w 0 (x) = e − x 2 4ε u 0 (x).Similarly, any solution of the SE or PWE can be transformed into a solution as well according to Interestingly, in connection with optical propagation, the modulation of the input function by e − x 2 4ε may signify lensing or Gaussian aperturing according to whether ε is real or purely imaginary [21].
Note that the modulating function in both Equations ( 27) and ( 28) is the fundamental solution Equations ( 16) or (18).Equation ( 28) confirms the Gaussian packets/beams as solutions to Equation ( 8) [3] and Equation ( 9) [4], following indeed from Equation ( 28) when setting ψ 0 (x) = 1, which amounts to ψ(x, ς) = 1 as well.In this case, the symmetry transformations simply manifest the symmetry of the solutions of Equation ( 10) with respect to the shift of the evolution variable ζ [21].In fact, Equations ( 27) and ( 28) follow directly from the additivity of the evolution operator U(ζ), being As said, the fundamental solution is the result of the evolution of the initial data δ(x); namely, referring to the HE, we can write Therefore, the Gaussian e − x 2 4ε can be regarded as resulting from δ(x) at the t = ε (apart from the factor 1 √ 4πε ), and hence This implies that In practice, it is as we had translated the origin of the "time" from t 0 = 0 to t 0 = −ε.
Of course, the same can be seen for Equation (28) with the simple replacement t → iς [21].
Central to the theory of the HE is the well-known Appell transformation [2,23].It establishes the possibility to pass from one solution of the HE to another according to the rule Apart from its role in connection with, for instance, the heat polynomials, the interest in the Appell transformation arises from its significative property, proven in [29], according to which it is essentially the only transformation mapping solutions of the HE into another in the sense that any simmetry transformation of the HE can be seen as composed by Appell transformations and scaling and shift of both coordinates x, t.
Appell transformation for the PWE (as well as for the SE) can be considered, of course, as discussed in [27,28], accordingly yielding the optical analogue of Equation (30) as Evidently, in conformity to the formal analogy of the equations of concern, the aforementioned result regarding the role of the Appell transformation in relation to the symmetry transformations of the HE [29] can be reformulated in connection with the optical Appell transformation.In fact, in [30] the latter has been proven to be the only symmetry transformation of the PWE in the sense that any symmetry transformation of the PWE can be obtained by composing Appell transformations with scaling and translations of both coordinates.
We conclude recalling that, as Equations ( 27) and ( 28) relate to the symmetry operator (for Equation ( 10)) e − 1 2ε K + = e − x 2 4ε , the Appell transformation, conveyed by Equations ( 30) and ( 31), relates to the symmetry operator represented by the Fourier trasform As said, the analysis here presented frames within an investigation aimed at establishing if, at which extent and in which form some properties of the evolution equations-in particular, those exemplified in connection with the just discussed HE, SE and PWE-might be extended to equations, which are not of evolution type or demand to deal with fractional differential operators, like the RHE, KGE,WE and RSE, mentioned in the Introduction.
As we will show in Sections 4 and 5, operator factorization à la Dirac can help in this.

Evolution Operator and Fractional Partial Differential Equations
As stressed in the previous discussion the evolution operator is a flexible tool which can be applied within a more general context going beyond the cases analyzed so far.For instance, it can effectively be applied to fractional differential equations, which, as is well known, are finding wide applications in modeling both physical and engineering systems [31][32][33][34][35].
Indeed, in the case of fractional evolution equations like the solution can be formally obtained as x .
The use of the Laplace transform identity [36] e where g α (ξ) is a Levy stable distributions or α-stable distributions, allows to write the solution of Equation ( 32) in the form [11] F On the other hand, the solution of the equation where O x is an operator acting on the x variable, can be obtained by the use of a slight redefinition of the evolution operator method.According to [37], where and γ(αr+1) is the Mittag-Leffler function.Furthermore the use of the identity [37] which is a fairly noticeable results, since it provides the solution of a fractional partial differential equation of the type of Equation (33) as the Laplace transform of its integer counterpart.If e.g., O x = ∂ 2 x and f (x) = e −x 2 , we find Finally the equation can be solved by merging the two procedures, thus finding It is evident that operational and integral transform methods are general and useful tools to deal with fractional derivative operators and also if fractional differential operators are involved as in the case of the relativistic heat equation, which reads whose solution is readily obtained in the form In this paper we will consider the problem of dealing with various forms of differential equations involving square and higher order root operators by the use of a generalization of the Dirac factorization method, which in our opinion may offer a new point of view to the theory of fractional operators.

Dirac-like Factorization to Disentangle Root Operator Functions
As is well known, the algebra of the operators is definitely different from that of c-numbers.Thus, for instance, the identity can not hold if A and B are numbers (real or complex).In contrast, it can hold if A and B are operators or matrices satisfying definite relations.We will refer to Equation (34) as Dirac-like factorization procedure.On reviewing already published results [12][13][14], we will firstly exemplify the procedure in the case of a square-root operator, and then will show how to extend the method to higher order root-operators.New results will then be established.

Square Root Function
The identity can hold if A and B are anticommuting operators or matrices, for which so This suggests to write down the square root function where α and β turn out to be "mathematical objects" such that in order for the desired equality, given by Equation ( 35), be satisfied.In principle, A and B can be either numbers or commuting operators; of course, the latter case is of major interest, and it will be considered here.
Clearly, α and β cannot simply be numbers; indeed, as a direct consequence of Equation ( 36), they must be traceless matrices with eigenvalues equal to ±1, and hence of order 2n × 2n, n = 1, 2, .., and determinant equal to (−1) n .Thus, on one side one looses the scalar nature of the original function, which in the stated identity Equation ( 35) is to be understood as a multiple of the 2n × 2n unit matrix I 2n×2n .In fact, Equation ( 35) conveys a matrix identity in a proper 2n-dimensional vector space, whose meaning and ultimate dimensions (following those of α and β) are dictated by the context inherent in the problem at hand.On the other side, one gains a root-free matrix form expression, that could facilitate the solution of the problem although it must be reintenpreted in the light of the gained degree (or, degrees) of freedom (naturally conveyed, as seen, by the procedure).In addition, the method can open new perspectives within the theory of fractional calculus, suggesting alternative formulations to already established treatments and/or definitions.
At a basic level, involving up to three addends in the square root, the smallest admissible dimension 2n = 2 is enough to ensure that the desired matrices α and β can be realized.So, on the basis of Equation ( 36), we may identify them with any two of the Pauli matrices We recall that I 2×2 being the 2 × 2 unit matrix, and ε jkl the Levi-Civita tensor: ) is an odd permutation of (1, 2, 3), 0 if any two indices are equal .

Eventually, we can write
or, more in general, As said, the correspondence of each of the operators involved in the square root with a specific Pauli matrix is a mere matter of convenience, possibly suggested by the problem under investigation.Therefore, the resulting matrix expression of the original (scalar) operator function is not unique.
Finally, we recall the well known property of the Pauli matrices that according to which the exponential matrix e aσ j belongs to the algebra spanned by the Pauli matrices and the unit matrix I 2×2 .

Estension of the Procedure to Higher-Degree Root Operator Functions
The obvious extension of the above procedure amounts to establishing whether it be possible to write down or, more in general, with m operators involved in, with the α and β matrices being identified through suitable conditions analogous to Equation (36).

Cube Root
Indeed, the factorization allowing for the disentanglement of the cube root operator as is possible for commuting operators A, B, and matrices α and β such to satisfy the three-term relations (in a sense paralleling the two-term relations in Equation ( 36)), which can equivalently be written as We see that α and β are traceless matrices, with eigenvalues conveyed by the third roots of unity: and hence Therefore, they must be of the order 3n × 3n, n = 1, 2, .., with determinant The matrices of smallest admissible dimension, provide a suitable pair of matrices satisfying the required conditions.They are seen to span, by repeated commutators, an 8-dimensional Lie algebra.For instance, For completeness' sake, we write down the expressions of the other τ-matrices, Each of them is such that τ 3 j = I 3×3 , and δ j = 1, j = 1, .., 8. Also, their commutators are the relevant structure constants f jkl being Of course, f jkl = −f kjl .Interestingly, only 4 commutators vanish, i.e., and accordingly, the involved pairs of matrices are the only ones that do not allow for the desired factorization à la Dirac of the sum of third-power operators.Note in fact that the relations in Equation ( 44) could not be satisfied by (non null) commuting matrices.Therefore, 24 possible pairs of matrices suitable for the disentanglement of the cube root are conveyed by the set of τ-matrices.
If a third term is added in the sum, amounting to the linearization issue a triplet of matrices is needed, such that each and each pair of them satisfy the relations in Equation ( 44), in addition to the further one the sum being over all the six possible products of the three matrices obtained from all their permutations p (∈ S 3 ).The relations in Equation (44) for each pair of matrices can as well be cast in the form as Equation ( 49) if applied to sets of three matrices in which two of them are the same.One can see that 24 suitable triplets of matrices can be extracted from the set of τ-matrices, the choice being a matter of convenience in conformity to the problem under analysis.
The τ-matrices have already been deduced in [43] in connection with the analysis of the fractional Dirac equation.There, the triplets of matrices allowing for Equation (48) are explicitly signalized.
We conclude by noticing that, as a consequence of that τ 3 j = I, the exponential matrix e aτ j turns out to be the sum of three terms; precisely, The coefficients A j (a), j = 0, 1, 2 are given by ; ; ; where 0 F 2 (•) denotes the generalized hypergeometric function.
Formally represented by the series [26,38] p F q (a 1 , .., a p ; b 1 , .., b q ; z) with (a) k ≡ γ(a + k)/γ(a) being the Pochhammer symbol, p F q converges for all finite z if p ≤ q.Hence, the A j 's are convergent series; also, they have been investigated in [39,40] as pseudo-hyperbolic functions.Equation ( 50) corresponds to Equation (40) on due account of the link between the hyperbolic and the hypergeometric functions: , ).
Finally, we note that the set {τ j } j=1,..,8 is closed with respect to the square, since the square of a τ-matrix just equals the matrix in the set with which it commutes; thereby, τ 2 1 = τ 8 , τ 2 2 = τ 7 , and so on.This property of the τ-matrices (i.e., the cube roots of the unit matrix) perfectly parallels that of the cube roots of unity; in fact, µ 2 0 = µ 0 , µ 2 1 = µ 2 , and As a consequence, the exponentials e aτ j belong to the algebra spanned (in general, over the complex field C) by the {τ j } j=1,..,8 and the unit matrix.

Quartic Root
Likewise, the disentanglement of the quartic root as demands for α and β matrices such to satisfy the four-term relations Therefore, we can say that the desired matrices α and β are traceless and with eigenvalues conveyed by the quartic roots of unity: As a consequence, they must be of the order 4n × 4n, n = 1, 2, .., with determinant The anticommutation relations in the second row of Equation ( 54) suggest a correspondence of the matrices α 2 , {α, β} and β 2 with σ-composed matrices.Thus, working with the smallest admissible dimension, i.e., 4n = 4, we start taking with √ σ 3 = 1 0 0 i .
By repeated commutators, we span a 15-dimensional Lie algebra of matrices ρ j j=1,..,15 such that ρ 4 j = I, ∀j.Let us explicitly write the other 12 matrices: Then, with the Lie bracket recast in the form we can specify the structure constants g jkl as with g jkl = −g kjl .However, not all the matrices in Equation (57) satisfy the condition on the determinant being in fact δ(ρ 5 ) = δ(ρ 8 ) = δ(ρ 11 ) = 1; accordingly, such matrices are not suitable for the desired factorization.Besides, they commute with each other, and each of them commute also with other four matrices in the set ρ j j=1,.., 15 .In turn, each of the other ρ-matrices commutes only with two elements in the set, one being ρ 5 , ρ 8 or ρ 11 .Of course, this can be deduced from the set of {j, k} indices that do not appear in the above reported table of structure constants g jkl .
An accurate analysis reveals that one can rely on 48 possible pairs of matrices allowing for Equation (53); as expected, the matrices ρ 5 , ρ 8 and ρ 11 are not inlcuded in any of the possible pairs.
Another property of the ρ-matrices is worth mentioning.It in a sense parallels that concerning the τ-matrices, and suggests a generalization to the m-th roots of the unit matrix.One can easily verify that the set ρ j j=1,..,15 ⊕ I 4×4 is closed under the square and the cube.Interestingly, we find that as a consequence, indeed, of that they are directly linked to the Pauli matrices.By applying the same terminolgy as for (complex) numbers, we can say that such 4 × 4 matrices are not primitive 4-th roots of the 4 × 4 unit matrix, being in fact also 2-th roots of the 4 × 4 unit matrix (In fact, µ 0 and µ 3 are not primitive 4-th roots of unity, whilst µ 1 and µ 2 are).This clearly clarifies why such matrices have determinant equal to 1.
In contrast, the squares of the other ρ-matrices just equal ρ 5 , ρ 8 or ρ 11 (apart from a minus sign), precisely that with which the matrix of concern commutes.Thus, for instance, ρ 2 1 = ρ 8 , ρ 2 2 = ρ 11 , ρ 2 3 = ρ 5 , and so on.Finally as to the cube, it is obvious that whereas the cube of any other ρ-matrix just equals the matrix (different from ρ 5 , ρ 8 and ρ 11 ) with which the matrix of concern commutes (apart from ±1 or ±i).Thus, for instance, ρ 3 1 = iρ 15 , ρ 3 2 = ρ 14 , ρ 3  3 = −ρ 13 , and so on.Then, as for the τ-matrices, we can say that the exponential matrix e aρ j belongs to the algebra spanned (in general, over the complex field C) by the ρ j j=1,..,15 and the unit matrix.In fact, the analog of Equation ( 50) can be written as with the coefficients B k (a), k = 0, 1, 2, 3 being correspondingly given by as an obvious extension of Equation ( 51) to the matrices Equation (57).Needless to say, for the non primitive 4-th root matrices ρ 5 , ρ 8 and ρ 11 we can write down an expression similar to Equation (40), namely e aρ l = cosh(a)I 4×4 + sinh(a)ρ l , l = 5, 8, 11.

m-th Roots
It is evident that with increasing m in Equation ( 41) the problem becomes even more complex.However, on the basis of the previous analysis, we can try to draw some basic issues, at least for the two-term case just displayed in Equation (41).
The identity in Equation ( 41) yields m + 1 relations involving terms of degree m in the mn × mn matrices α and β.Firstly, the latter come to be the m-th roots of the unit matrix, being and hence their eigenvalues can be written as It is easy to see that we can say that the determinant of the matrices α and β is 1) .
As to the other conditions ensuring Equation (41) be satisfied, they can be synthesised in the form p,q where the sum is intended to comprise all the powers p j and q i such that j p j = l and i q i = k, for any choice of integers (l, k) such that l = 0, k = 0 and l + k = m, thus yielding m! l!k! terms of degree m.
On the basis of the previous analysis, we can also state that in general the set {ν} j=1,..,m 2 −1 ⊕I m×m is closed under the powers 2, 3, ...m.Consequently, we can write the group element e aν j in the general form where the coefficients

Possible Applications: Practical and Conceptual Issues
The Dirac-like factorization procedure can be applied to various (physical and/or mathematical) contexts, and also be variously finalized.
We will illustrate some applications in connection with root functions of differential operators, specifically, square and cube roots respectively of second-order and third-order differential operators.

Square Root of Differential Operators
As said, the factorization approach to square root of differential operators can conveniently be exploited in connection with evolution equations ruled by pseudo-differential operators [12,14] as well as within the theory of fractional calculus, where it may open new perspectives [13,14].
Let us discuss both issues in some detail.
We will refer to it as relativistic-like free evolution equation.In fact, setting ζ = i ct λc and ξ = x λc , λ c = h mc being the Compton wavelength of the particle, the above equation would yield the (1+1)D version of the Salpeter Equation ( 5), and accordingly ψ(ξ, ζ) would represent the particle wave function.However, we will limit to consider in more detail real evolution variables ζ, so that the HE naturally arises as the "nonrelativistic" counterpart of Equation (62).A few comments on the solution à la Dirac of the (1+1)D Salpeter equation will be however given later.The Dirac-like "linearization" procedure turns the problem of the solution of Equation ( 62) into that of the solution of the system of two coupled linear homogeneous first order partial differential equations The solution to the above "evolution equation" is indeed immediately written in the form in full analogy with Equation (12).Whatever be the specific Pauli matrices chosen in the factorization, on account of that σ 2 j = I 2×2 and {σ j , σ k } = 2δ j,k I 2×2 , the direct evaluation of the exponential function in Equation (64) [42] would yield for the evolution matrix the explicit expression Let us apply such a result to specific initial data.Indeed, we consider the input vector Evidently, the definite expression of the evolved vector depends on the specific choice of the Pauli matrices in the square-root factorization.
Thus, with the evolution of the vector ψ is dictated by C and S being the cosh-and sinch-operator functions entering Equation (65), i.e., Therefore, with ψ 0 (ξ) given by Equation (66), we obtain Equivalently, resorting to the Fourier transform, defined as the vector ψ(ξ, ζ) would be given by The functions C and S represent the Fourier images of the operators C and S , namely Equations ( 68) and (70) yield the same results for the two vector components, which are plotted in Figure 1 vs. ξ at some ζ.   3 shows the behavior of the squared moduli of the two vector components ψ 1 and ψ 2 , obtained as Dirac approach based solutions of the (1+1)D Salpeter equation for the input given by Equation (66).As said, the latter follows from Equation (62) with the replacement ζ → iτ, τ = ct λc .Of course, this amounts to replacing the hyperbolic functions with the circular functions, i.e., cosh → cos and sinh → i sin, in all the relevant expressions.Accordingly, complex expressions for ψ 1 and ψ 2 are obtained, thus demanding to plot the respective squared moduli.Finally, as in Figure 2, Figure 4a shows the Euclidean norm ||ψ| | 2 = ψ 2 1 + ψ 2 2 of the vector solution ψ at the τ's pertaining to Figure 3.It can be compared with the squared amplitude of the wave function obtained as solution of the Salpeter equation through the Fourier transform based approach for the Gaussian input Equation (66), namely ψ 0 (ξ) = e − ξ 2 4 .Such an approach yields the solution of the (1+1)D Salpeter equation in the form [10,11] corresponding to the initial wave function ψ 0 (ξ), whose Fourier transform is ψ 0 (κ).

Suggesting Alternative Formulations in Fractional Calculus
Another possible context of application of the Dirac-like "linearization" procedure is that of the theory of the fractional calculus [43][44][45].As an example, let us consider the operator a being an arbitrary constant.For it the following interpretation, can be worked out, resorting to the integral representation of the operator L −ν , (ν) > 0, as which reproduces for operators the well-known Laplace-transform identity for c-numbers.Note that the shift operator e −s∂x under the integral yields e −s∂x f (x) = f (x − s).
Alternatively, in the light of the above analysis, the operator can be replaced by the operator matrix for any specific choice of the inherent Pauli matrices, thus opening new perspectives within the theory of fractional calculus.The operator nature of the l.h.s.would be conveyed by the matrix nature of r.h.s.; indeed, √ ∂ x may be seen as acting on 1, thus giving according to the Euler definition of fractional derivative x µ−ν .
Thus, the operator in Equation ( 72) can be regarded as acting in a 2D vector space through the matrix or also through , each matrix being obtained in correspondence with a specific choice of the Pauli matrices in Equation (74).In our opinion, the view conveyed by the above analysis deserves to be explored.

Cube Root of Differential Operators
Paralleling the analysis, developed in connection with the square root functions of operators, regarding the evolution Equation (62), we may consider an evolution equation involving a cube root of the differential operators ∂ 3 ξ , namely ψ(ξ, 0) = ψ 0 (ξ).
By the Dirac-like procedure, it is recast into the system of three coupled linear homogeneous partial differential equation of the first order for the three component vector ψ(ξ, 0), for any choice of the suitable pairs of the afore-introduced τ-matrices.
As before, the solution can be formally written as However, in order to get an explicit expression for the evolution matrix one needs to resort to appropriate ordering techniques since the τ-matrices involved in the linearization of the original cube root operators should not commute with each other.Just to fix mind, let us work with the matrices in Equation (46).Then, we apply the Zassenhaus formula [46,47] giving the exponential of the sum of two operators as the in general infinite product of operators according to where the first terms in the product explicitly write as Then, in the case of Equation ( 78) we find that thus enabling us to write for U (ξ, τ) at the third order in the evolution parameter ζ the expression the latter expression being allowed by the commutation [τ 4 , τ 5 ] = 0.In turn, each exponential can be written in the form as Equation (50).

Relativistic-like Evolution Equation: A "Direct" Solving Method
We will reconsider Equation (62), rewritten here for convenience's sake, ψ(ξ, 0) = ψ 0 (ξ), and approach its solution through a more "direct" method.Equation (62) has been the object of the recent analysis in [48].Here, we will further elaborate that analysis, trying to reproduce as much as possible the solving procedure inherent to the HE.Interestingly, we will find many analogies.

Evolution Operator
In full analogy with Equation (10), the solution to Equation (79) can be expressed in terms of the evolution operator, i.e., with U(ζ) being explicitly given by the exponential operator Evidently, it reproduces the evolution operator in Equation ( 14) pertaining to Equation (10) to the lowest order in ∂ 2 ξ of the power series expansion of 1 − ∂ 2 ξ , apart from some scaling and multiplying factors, as we will detail later.

Evolution Operator as McDonald Transform
In order to determine the action of U(ζ) on the initial data ψ 0 (ξ), we resort to the result [49][50][51] setting in particular a = 1 and b One can recognize under the integral the "non-relativistic" evolution operator in Equation ( 14) (at "time" ζ 2 4s 2 ).Accordingly, on account of Equation ( 15), we can write It is convenient to change the order of the integrals, and then to exploit the integral representation of the modified Bessel function of the second kind (or, McDonald function) K ν , i.e., [26,51,52] Thereby, we end up with the integral transform expression for the solution of Equation (79): Clearly, the same result can be obtained working in the Fourier domain.In fact, since 1 − ∂ 2 ξ e ±iκξ = √ 1 + κ 2 e ±iκξ , Equation (79) signifies in the Fourier conjugate κ-domain.Here, ψ(κ, ζ) denotes the Fourier transform (see Equation ( 69)) of ψ(ξ, ζ) with respect to ξ.
The solution to Equation ( 87) is then easily formulated as which, when transformed back to the ξ-domain, on account of yields Equation (86) for ψ(ξ, ζ).This in turn legitimates the assumption underlying Equation (83).Equation (86) conveys the interesting result that, just as for Equation (10) governed by the non-relativistic free-Hamiltonian, written in Equation (11), the solution to Equation ( 79) is given by an integral transform of the initial data, having the form of a convolution product, whose kernel involves the McDonald function K 1 (z).Accordingly, we will refer to Equation (86) as McDonald transform.As basic properties of the McDonald function K ν (z), we may recall that it is real when ν is real and z is positive, and that K ν (z) = K −ν (z).
Further analogies between Equations ( 10) and (79) can be drawn.Firstly, in the light of Equation (88), we see that thus enabling us to refer to V (ξ, ζ) as the fundamental solution of Equation (79).Correspondingly, the opposite limit ζ 1 and |ξ − ξ | ζ yields the "non-relativistic" solution.In fact, we can exploit the asymptotic properties of K 1 , giving [26,52] according to which we get Then, we replace the square root √ ζ 2 + ξ 2 = ζ 1 + ξ 2 ζ 2 by its power series expansion in ξ 2 ζ 2 at zero order in the denominator, thus yielding ζ, and at the first order in the exponent, giving ζ+ ξ 2 2ζ .Eventually, we end up with the expression which reproduces the Poisson transform in Equation ( 15) up to the term e −ζ .The latter accounts in a sense for the "rest energy", and it does not appear in Equation (15).Also, the evolution variable is half that pertaining to Equation ( 10) since the non-relativistic limit of Equation (79) would yield 1 2 ∂ 2 ∂ξ 2 ; the factor 1  2 is inglobed in the variable ζ in Equation ( 10), emerging indeed in the specific cases, for instance, of the SE and the PWE.

Symmetry Transformations and "Polynomial" Solutions
A symmetry property analogous to Equation ( 27) for the case u(x, t) = 1, in practice the corresponding relation of the evolution of the fundamental solution for the case of concern can be established in the case of concern as well.
In fact, on the basis of the relation for the fundamental solution Equation (89), we can immediately verify that with we obtain explicitly yielding the solution As remarked in Section 2.3, the above simply manifests the symmetry of the solutions of Equation ( 79) with respect to a shift of the evolution variable ζ.In both cases, the fundamental solution evolves through a simple shift of the evolution variable.
The possibility of obtaining an expression corresponding to Equation ( 27) in its general form is under investigation.
Figure 5 compares the evolution of the Gaussian ψ 0 (ξ) = e − ξ 2 4 (in practice, the fundamental solution of the HE) as it should occur in accord with Equation (10) 4(1+ζ) , and in accord with the relativistic Equation (79), i.e., through Equation (86).The (ζ, ξ)-contourplots are shown in both cases; in the latter case, as a check, the wavefunction at ζ > 0 has also been obtained through the Fourier transform method.
The "relativistic" evolution shows a quicker attenuation of the wavefunction amplitude than the "non-relativistic" one.In order to have a function, which, just as the Gaussian e − ξ 2 4 , takes on the unit value at ξ = 0, the factor π K 1 (ε) is added to Equation (93), so that the function is plotted in Figure 6.
For completeness'sake, Figure 7 plots the profiles of the wavefunctions of concern in the same panel at different ζs.We see that, with increasing ε, ψ D (ξ, ζ, ε) tends to be even more similar to the ψ G (ξ, ζ)| nr .Notably, as corresponding to initial monomials can be found for the HE, SE and PWE as recalled in Section 2.3, explicit solutions corresponding to initial monomials can be found for Equation (79) as well.
Furthermore, as the heat polynomials realize the power series expansion of the evolved form (through Equation ( 7)) of the eigenstates e χξ of the derivative operator ∂ ∂ξ , i.e., e χξ e χ 2 ζ , the relativistic heat polynomials realize the power series expansion of the corresponding "relativistically" evolved form, i.e., for |χ| ≤ 1.
The above relies on the Neumann-type series for the MacDonald function (also referred to as multiplication theorem) [26,52]  Therefore, the generating functions pertaining to the p n s as well as the v n s can be ascribed to the Appèl family [53,54].Also, according to [54], we can express the relativistic heat polynomials in the operator form Since the operator 1 − ∂ 2 ξ −1/2 can be understood as one recovers the relations in Equation ( 23) for the v n s as well as the inherent differential Equation ( 24), taking only the 0-th order term in the series.Also, we see that the raising operator for the p n s is a ξ-differential operator of order 2 n−1 2 + 1, and accordingly Equation ( 98) is of order 2 n 2 .Finally, it is worth noting that the p n s belong to the family of non-local polynomials, since they obey the integro-differential Equation (98) instead of an ordinary differential equation like the v n s.In fact, on the basis of the representations in Equations ( 73 Other properties of Equation (86) are under investigation.

Conclusions
As a particular aspect of the general issue of the analysis of equations involving fractional or pseudo-differential operators, we have reviewed the square-root operator factorization method à la Dirac along with its extension to higher-degree root operators, as recently suggested in the literature [12][13][14].A deeper analysis of the cube and quartic root operators has been presented along with a precise characterization of the Lie algebras of the pertinent matrices.
In addition, evolution equations ruled by square root operator functions, referred to as relativistic-like free evolution equations, have been considered, further elaborating the analysis developed in [48].In fact, a closed form expression of the solution has been deduced as an integral transform of the initial data.
We have referred to it as McDonald transform since the kernel involves the McDonald function K 1 .The presentation here has been finalized to a comparison between the properties of the non-relativistic and relativistic equations, in order to establish which properties of the former can be recovered to the latter.In particular, fundamental solution with its evolution property and solutions arising from initial monomials have been found for the relativistic evolution equation in full analogy with the non relativistic one.Further investigations of the McDonald transform have been suggested.

Figure 2
Figure 2 shows the Euclidean norm |ψ| 2 = ψ 2 1 + ψ 2 2 of the vector ψ at the ζ's pertaining to Figure 1.Paralleling Figure 1, Figure3shows the behavior of the squared moduli of the two vector components ψ 1 and ψ 2 , obtained as Dirac approach based solutions of the (1+1)D Salpeter equation for the input given by Equation (66).As said, the latter follows from Equation (62) with the replacement ζ → iτ, τ = ct λc .Of course, this amounts to replacing the hyperbolic functions with the circular functions, i.e., cosh → cos and sinh → i sin, in all the relevant expressions.Accordingly, complex expressions for ψ 1 and ψ 2 are obtained, thus demanding to plot the respective squared moduli.
calculus.Note that the l.h.s.can be interpreted as the Laplace transform of f (s) = e − b s / √ s.Then, assuming that the identity expressed by Equation (82) holds also when the parameters a, or b or both are (commuting) operators, we can write