# Remarks on Root Operators and Extensions of the Dirac Factorization Method

Julius Lehmann

► To cite this version:

Julius Lehmann. Remarks on Root Operators and Extensions of the Dirac Factorization Method. 2026. hal-05470467

HAL Id: hal-05470467

https://hal.science/hal-05470467v1

Preprint submitted on 21 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Remarks on Root Operators and Extensions of the Dirac Factorization Method

Julius Lehmann*

*Physics of Complex Biosystems, Technical University of Munich, 85748 Garching, Germany*

January 21, 2026

## Abstract

In this work, we explore generalized Dirac-inspired factorizations of a variety of differential equations commonly encountered in mathematical physics, highlighting how different choices of root operators reveal hidden algebraic structures and physical interpretations. Throughout, we emphasize the flexibility of the factorization approach, illustrating how it allows for the construction of alternative formulations that may be fully linear, retain time or spatial memory, or circumvent the explicit appearance of fractional derivatives. We also provide a general method for determining coefficients in Clifford algebra expansions. These results highlight the utility of algebraic factorization techniques for uncovering hidden structures, clarifying phenomenology, deepening our understanding of their underlying dynamics and offering novel perspectives on classical and fractional differential equations in mathematical physics.

**Keywords:** Dirac factorization; Fractional differential operator; Ordinary and partial differential equation; Pseudo-differential operators

## Introduction

Already at the advent of calculus, factorization methods proved pivotal in providing solutions to higher-order differential equations by decomposing them into first-order differential operators. One such example is d'Alembert's solution to the one-dimensional wave equation (also derived by him) [5]:

$$ \left( \frac{\partial^2}{\partial x^2} - \frac{1}{c^2} \frac{\partial^2}{\partial t^2} \right) \phi(x, t) = 0, $$

which he obtained by factoring the equation into two components using the interchangeability of partial derivatives:

$$ \left( \frac{\partial}{\partial x} - \frac{1}{c} \frac{\partial}{\partial t} \right) \left( \frac{\partial}{\partial x} + \frac{1}{c} \frac{\partial}{\partial t} \right) \phi(x, t) = 0, $$

<footnote>*e-mail: julius.lehmann(at)tum.de</footnote>
from which we can readily read off the solution as the sum of two contributions,

$$ \phi(x, t) = \phi_1(x - ct) + \phi_2(x + ct), $$

namely a right-traveling function and a left-traveling function, which can be further specified through the imposition of initial conditions. The next major impetus occurred in the 1920s, when Dirac employed his now eponymous factorization method to lay the foundations of relativistic quantum mechanics, while simultaneously providing a rigorous justification of electron spin and the theoretical prediction of antimatter [10, 9, 14]. From that point onward to the present day, numerous works have investigated how the Dirac factorization method may reveal underlying (symmetry) structures and hidden physics, particularly when addressing pseudo-differential operators and within the framework of fractional calculus [2, 7, 6, 16, 11, 20]. Notably, this approach yields distinct perspectives on various “classical” equations in mathematical physics, offers novel physical interpretations, and demonstrates how specific phenomenologies emerge naturally within this framework [6, 11].

In the present work, we explore the application of the Dirac factorization method to partial and ordinary differential equations, possibly involving fractional or pseudo-differential operators. We demonstrate that generalizing the conventional approach to more elaborate “roots”, that is, factors which do not simply correspond to the rooted components of the differential operator but instead incorporate higher- or lower-order operators and combinations thereof, gives rise to new insights and phenomenologies. We consider a range of examples, beginning with the wave equation and demonstrating how this methodology leads to a one-way wave equation. We then explore how different choices of roots for the telegrapher equation give rise to, in part, diametrically different physical interpretations, and further examine how fractional-order differential equations can be transformed into higher-order equations and how new physical insights can be extracted from this procedure.

### Example 1: Wave Equation

We begin our analysis with the historically most significant case, namely the wave equation. In three dimensions, it can be written as the second-order partial differential equation

$$ \left( \Delta - \frac{1}{c^2} \frac{\partial^2}{\partial t^2} \right) \phi(x, t) = 0, \qquad (1) $$

where the differential operator is also known as the d’Alembertian. Depending on the physical interpretation, the function $\phi(x, t)$ may represent different physical quantities, such as mechanical waves (e.g., water waves, sound waves, or seismic waves) or electromagnetic waves (including light waves). To introduce the core idea underlying the Dirac factorization method, we consider the simple example of the one-dimensional wave equation

$$ \left( \frac{\partial^2}{\partial x^2} - \frac{1}{c^2} \frac{\partial^2}{\partial t^2} \right) \phi(x, t) = 0, \qquad (2) $$

for which we seek to compute the corresponding root operator

$$ \sqrt{\frac{\partial^2}{\partial x^2} - \frac{1}{c^2} \frac{\partial^2}{\partial t^2}} \phi(x, t) = 0. \qquad (3) $$
In principle, one could treat this pseudo-differential operator through its power series expansion or by employing integral transform techniques; however, such approaches generally lead to non-local or physically intractable equations. Instead, following Dirac's original reasoning, we suppose that the root can be expressed as a linear combination of linear partial derivatives, namely

$$ \sqrt{\frac{\partial^2}{\partial x^2} - \frac{1}{c^2} \frac{\partial^2}{\partial t^2}} \phi(x, t) = \left(A \frac{\partial}{\partial x} + \frac{i}{c} B \frac{\partial}{\partial t}\right) \phi(x, t), \quad (4) $$

for some unknown coefficients A and B. By squaring both sides, we obtain a set of algebraic conditions that these coefficients must satisfy,

$$ A^2 = 1, \quad B^2 = 1, \quad AB + BA = 0. \quad (5) $$

Because scalar-valued coefficients always commute, no such numbers can satisfy these relations. To circumvent this obstacle, we instead consider higher-order quantities, namely matrices. The minimal set of matrices that fulfills the conditions in Eq. (5) is provided by the three Pauli matrices [18],

$$ \sigma_1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad \sigma_2 = \begin{pmatrix} 0 & i \\ -i & 0 \end{pmatrix}, \quad \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}. \quad (6) $$

Using these matrices, we can write the root of the one-dimensional wave operator in the form

$$ \left( \sigma_i \frac{\partial}{\partial x} + \frac{i}{c} \sigma_j \frac{\partial}{\partial t} \right) \vec{\phi}(x, t) = \vec{0}, \quad (i \neq j) \quad (7) $$

where the price paid for eliminating the pseudo-differential operator is the introduction of a higher-dimensional wave function. Since there are six possible ordered combinations of the indices (1, 2, 3), the resulting equations are not unique and, in fact, depend explicitly on the particular choice of indices i and j used to realize the factorization. The complete list of distinct equations corresponding to all admissible pairs (i, j) is given by

$$ \begin{array}{l} (1,2): \left\{ \begin{array}{l} \displaystyle \frac{\partial \phi_2}{\partial x} + \frac{1}{c} \frac{\partial \phi_2}{\partial t} = 0 \\ \displaystyle \frac{\partial \phi_1}{\partial x} - \frac{1}{c} \frac{\partial \phi_1}{\partial t} = 0, \end{array} \right. \qquad (1,3): \left\{ \begin{array}{l} \displaystyle \frac{\partial \phi_2}{\partial x} + \frac{i}{c} \frac{\partial \phi_1}{\partial t} = 0 \\ \displaystyle \frac{\partial \phi_1}{\partial x} - \frac{i}{c} \frac{\partial \phi_2}{\partial t} = 0, \end{array} \right. \qquad (2,3): \left\{ \begin{array}{l} \displaystyle \frac{\partial \phi_2}{\partial x} - \frac{1}{c} \frac{\partial \phi_1}{\partial t} = 0 \\ \displaystyle \frac{\partial \phi_1}{\partial x} - \frac{1}{c} \frac{\partial \phi_2}{\partial t} = 0. \end{array} \right. \quad (8) \end{array} $$

Each of these equations gives rise to a distinct phenomenology. In particular, the case $(i, j) = (1, 2)$ reproduces the original d'Alembert factorization into forward- and backward-propagating wavelets, whereas the remaining combinations exhibit more complex dynamical behavior. We illustrate the solutions for these unique index choices in Figure 1 for a Gaussian initial condition.

This procedure can, of course, be extended to higher-dimensional wave equations. In two dimensions, we consider the partial differential equation

$$ \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} - \frac{1}{c^2} \frac{\partial^2}{\partial t^2} \right) \phi(x, y, t) = 0, \quad (9) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Solutions to the different root operators of the one-dimensional wave equation (8) for a Gaussian initial condition $\phi_1(x, 0) = e^{-x^2}$, $\phi_2(x, 0) = 0$. The dashed black Gaussian represents the initial state, while progressively darker color hues indicate states at later times. Different choices of $(i, j)$ lead to distinct phenomenologies, with (1, 2) corresponding to the splitting into forward- and backward-propagating waves, and (2, 3)/(1, 3) yielding a (complexified) destructive and constructive wave front.</figcaption></figure>

for which the corresponding root operator takes the form

$$ \left( \sigma_i \frac{\partial}{\partial x} + \sigma_j \frac{\partial}{\partial y} + \frac{i}{c} \sigma_k \frac{\partial}{\partial t} \right) \vec{\phi}(x, y, t) = \vec{0}, \quad (i \neq j \neq k) \quad (10) $$

thereby employing the full set of Pauli matrices. In the historically most relevant case of the full three-dimensional wave equation,

$$ \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} - \frac{1}{c^2} \frac{\partial^2}{\partial t^2} \right) \phi(x, y, z, t) = 0, \quad (11) $$

we find the associated root operator to be

$$ \left( \gamma_i \frac{\partial}{\partial x} + \gamma_j \frac{\partial}{\partial y} + \gamma_k \frac{\partial}{\partial z} + \frac{i}{c} \gamma_\ell \frac{\partial}{\partial t} \right) \vec{\phi}(x, y, z, t) = \vec{0} \quad (i \neq j \neq k \neq \ell), \quad (12) $$

which is more commonly recognized as the Dirac equation for relativistic spin-$\frac{1}{2}$ particles. The matrices introduced here are the gamma, or Dirac, matrices, which provide a matrix
<figure><img src="image_3.png" /><figcaption>Figure 2: Solution to the two-dimensional one-way wave equation for an initial Gaussian state $\phi_1(x, y, 0) = e^{-x^2-y^2}$, $\phi_2(x, y, 0) = 0$, with the top row displaying the full two-dimensional height field of $\phi_1$ and the bottom row showing a slice along the x-axis. Superimposed on the slice is the solution to the standard two-dimensional wave equation, plotted as a black dotted line. We observe that the dominant wavefront propagates exclusively in the +x-direction.</figcaption></figure>

representation of the Clifford algebra $Cl_{4,0}(\mathbb{C})$. The emergence of the Clifford algebra in this context is not accidental but is instead deeply connected to rotation groups, though we shall not elaborate on that connection here.

Considering the two-dimensional case, we can isolate the x-derivative by multiplying with $\sigma_i$, since this matrix is involutory:

$$ \frac{\partial \vec{\phi}}{\partial x} + \left( \sigma_i \sigma_j \frac{\partial}{\partial y} + \frac{i}{c} \sigma_i \sigma_k \frac{\partial}{\partial t} \right) \vec{\phi}(x, y, t) = \vec{0}. \qquad (13) $$

Because the product of two Pauli matrices is itself proportional to a Pauli matrix, we find that this equation describes waves with wavenumbers ($k_x = \pm\sqrt{\omega^2 - k_y^2}, k_y$) and frequencies $\omega$ for suitable choices of the indices $(i, j, k)$. We can therefore interpret this relation as the dispersion relation of two-dimensional waves propagating exclusively in the +x- or -x-direction, and hence as a two-dimensional one-way wave equation. A convenient choice of indices is $(i, j, k) = (1, 3, 2)$, which yields the simple system of equations

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \phi_1}{\partial x} - \frac{1}{c} \frac{\partial \phi_1}{\partial t} - \frac{\partial \phi_2}{\partial y} = 0 \\[1.5ex] \displaystyle \frac{\partial \phi_1}{\partial y} + \frac{\partial \phi_2}{\partial x} + \frac{1}{c} \frac{\partial \phi_2}{\partial t} = 0. \end{array} \right. \qquad (14) $$

This system describes waves propagating in the +x-direction for initial conditions of the form $\vec{\phi}(x, y, 0) = (\phi_0(x, y), 0)$, with propagation carried by the component $\phi_1$. In Fig. 2, we illustrate the evolution of an initial Gaussian peak governed by Eq. (14).

Thus far, we have factored the second-order wave equation into systems of linear equations in space and time. We now demonstrate how this approach generalizes to more elaborate roots. The Dirac factorization method allows us to decompose the wave equation into different lower-order equations that may combine higher- and lower-order
derivatives. To illustrate such a construction, we consider a “diffusion”-type factorization of the wave equation, in which we retain second-order spatial derivatives while having the time derivative in first order. For simplicity, we restrict the discussion to the two-dimensional wave equation, although the procedure extends straightforwardly to higher-dimensional cases. We begin once more with the standard two-dimensional wave equation,

$$ \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} - \frac{1}{c^2} \frac{\partial^2}{\partial t^2} \right) \phi(x, y, t) = 0, \qquad (15) $$

and introduce the following ansatz for a possible factorization:

$$ \left( \omega_1 + \omega_2 \frac{\partial^2}{\partial x^2} + \omega_3 \frac{\partial^2}{\partial y^2} + \frac{i}{c} \omega_4 \frac{\partial}{\partial t} \right) \vec{\phi}(x, y, t) = \vec{0}. \qquad (16) $$

Requiring that the square of this operator reproduces the original wave equation yields the following set of algebraic conditions for the matrices $\omega_i$:

$$ \begin{array}{l} \omega_4^2 = 1, \quad \omega_2^2 = \omega_3^2 = \omega_1^2 = 0, \quad \{\omega_4, \omega_1\} = \{\omega_4, \omega_2\} = \{\omega_4, \omega_3\} = \{\omega_2, \omega_3\} = 0 \\ \{\omega_2, \omega_1\} = \{\omega_3, \omega_1\} = 1. \end{array} \qquad (17) $$

We can satisfy this set of conditions by constructing the matrices as linear combinations of Pauli or Dirac matrices:

$$ \omega_4 = \gamma_i, \quad \omega_2 = \omega_3 = \frac{\gamma_j + i\gamma_k}{2}, \quad \omega_1 = \frac{\gamma_j - i\gamma_k}{2}, \qquad (18) $$

where we employed the algorithm outlined in Appendix B. By using the identity $\omega_4^{-1} = \omega_4$, we can rewrite this “diffusion” factorization in a more compact form by isolating the time derivative. In this way, we transform the wave equation into a diffusion equation with tensorial, and hence anisotropic, diffusion coefficients as well as a non-isotropic source term:

$$ \frac{\partial \vec{\phi}}{\partial t} = i c \omega_4 \left( \omega_1 + \omega_2 \frac{\partial^2}{\partial x^2} + \omega_3 \frac{\partial^2}{\partial y^2} \right) \vec{\phi}(x, y, t). \qquad (19) $$

Although we do not further analyze the precise physical implications of this equation, it suggests that wave phenomena may equivalently arise from diffusive dynamics of the components of a state vector under appropriate initial conditions. Moreover, for specific choices of gamma matrices and by restricting the components to be strictly real, we obtain, upon focusing on the first two components,

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \phi_1}{\partial t} = -\frac{c}{2} \Delta \phi_2 + \frac{c}{2} \phi_2 \\ \\ \displaystyle \frac{\partial \phi_2}{\partial t} = \frac{c}{2} \Delta \phi_1 - \frac{c}{2} \phi_1. \end{array} \right. \qquad (20) $$

By introducing the complex-valued wave function $\psi = \phi_1 + i\phi_2$, we can combine these equations and find that the wave equation factorizes into a set of Schrödinger equations. When carried further, this observation suggests that the two-dimensional Schrödinger equation may be interpreted as a “diffusion”-type factorization of an underlying wave equation, conceptually consistent with the foundations of quantum mechanics, namely

$$ \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} - \left( \frac{2m}{\hbar} \right)^2 \frac{\partial^2}{\partial t^2} \right) \phi(x, y, t) = 0, \qquad (21) $$
with an effective “propagation speed” given by $c = \frac{\hbar}{2m}$. We emphasize, however, that the phenomenological implications of this equation differ from those of the Schrödinger equation.

Finally, we present two additional factorizations without further discussion. The first combines transport behavior in the x-component with diffusive behavior in the y-component:

$$ (\omega_1 + \omega_2 \frac{\partial}{\partial x} + \omega_3 \frac{\partial^2}{\partial y^2} + \frac{i}{c} \omega_4 \frac{\partial}{\partial t}) \vec{\phi}(x, y, t) = \vec{0}, \qquad (22) $$

with one particular solution given by

$$ \omega_4 = \gamma_i, \quad \omega_2 = \gamma_j, \quad \omega_3 = \frac{\gamma_k - i\gamma_\ell}{2}, \quad \omega_1 = \frac{\gamma_k + i\gamma_\ell}{2}. \qquad (23) $$

The second retains wave-like behavior in the y-component while introducing friction-like terms through the x-dependence:

$$ (\omega_1 + \omega_2 \frac{\partial}{\partial x} + \omega_3 \frac{\partial^2}{\partial y^2} - \frac{1}{c^2} \omega_4 \frac{\partial^2}{\partial t^2}) \vec{\phi}(x, y, t) = \vec{0}, \qquad (24) $$

with one possible realization given by

$$ \omega_4 = \frac{\gamma_i - i\gamma_j}{2}, \quad \omega_2 = \gamma_k, \quad \omega_3 = \frac{\gamma_i - i\gamma_j}{2}, \quad \omega_1 = \frac{\gamma_i + i\gamma_j}{2}. \qquad (25) $$

In light of these constructions, one might ask whether it is possible to factor the wave equation into two completely separate parts, one governing the x-dependence and the other governing the y-dependence, namely

$$ (\omega_1^{(x)} + \omega_2^{(x)} \frac{\partial}{\partial x} + \frac{i}{c} \omega_3^{(x)} \frac{\partial}{\partial t})^2 (\omega_1^{(y)} + \omega_2^{(y)} \frac{\partial}{\partial y} + \frac{i}{c} \omega_3^{(y)} \frac{\partial}{\partial t})^2 \vec{\phi}(x, y, t) = \vec{0}. \qquad (26) $$

However, such a factorization leads to a contradiction, as the matrices would need to be simultaneously nilpotent and invertible. We find that this contradiction persists for all attempts to separate the x- and y-dependence across the combinations we have explored.

### Example 2: Diffusion Equation

The diffusion equation with diffusion constant $\alpha$,

$$ \left(\frac{\partial}{\partial t} - \alpha \frac{\partial^2}{\partial x^2}\right) \phi(x, t) = 0, \qquad (27) $$

provides useful pedagogical insight into how different factorizations and choices of state-vector components give rise to distinct physical interpretations. In one spatial dimension, we already know that its linearization yields the standard continuity equation together with Fick's law,

$$ \left\{ \begin{array}{l} \frac{\partial \phi}{\partial t} + \frac{\partial j}{\partial x} = 0 \\ j = -\alpha \frac{\partial \phi}{\partial x}. \end{array} \right. \qquad (28) $$
To verify this intuition and to demonstrate that these relations indeed follow from the factorization procedure, we write the root operator in the form

$$ \left(\frac{\partial}{\partial t} - \alpha \frac{\partial^2}{\partial x^2}\right) \vec{\phi}(x,t) = \left(\eta_0 \frac{\partial}{\partial t} + i\sqrt{\alpha} \eta_1 \frac{\partial}{\partial x} + \eta_c\right)^2 \vec{\phi}(x,t) = \vec{0}, \quad (29) $$

from which the following algebraic relations immediately follow:

$$ \eta_0^2 = 0, \quad \eta_1^2 = 1, \quad \eta_c^2 = 0, \quad \{\eta_0, \eta_1\} = 0, \quad \{\eta_0, \eta_c\} = 1, \quad \{\eta_1, \eta_c\} = 0. \quad (30) $$

A particular realization of these relations is found using the algorithm outlined in Appendix B and again provided by linear combinations of the Pauli matrices,

$$ \eta_0 = \frac{\sigma_i \pm i\sigma_j}{2}, \quad \eta_1 = \sigma_k, \quad \eta_c = \frac{\sigma_i \mp i\sigma_j}{2}. \quad (31) $$

Setting $(i, j, k) = (2, 1, 3)$, we obtain the pair of equations

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \phi_2}{\partial t} - \sqrt{\alpha} \frac{\partial \phi_1}{\partial x} = 0 \\[1em] \displaystyle \sqrt{\alpha} \frac{\partial \phi_2}{\partial x} - \phi_1 = 0, \end{array} \right. \quad (32) $$

which, upon identifying $\phi_2 \to \phi$ and $\phi_1 \to -j/\sqrt{\alpha}$, reproduces the standard continuity equation and Fick's law. While this interpretation is straightforward, it naturally raises the question of which physical interpretations arise from other choices of the index triple $(i, j, k)$. For instance, selecting $(i, j, k) = (2, 3, 1)$ yields

$$ \left\{ \begin{array}{l} \displaystyle \phi_1 - \phi_2 - \frac{\partial \phi_1}{\partial t} - \frac{\partial \phi_2}{\partial t} + 2\sqrt{\alpha} \frac{\partial \phi_2}{\partial x} = 0 \\[1em] \displaystyle \phi_1 - \phi_2 + \frac{\partial \phi_1}{\partial t} + \frac{\partial \phi_2}{\partial t} + 2\sqrt{\alpha} \frac{\partial \phi_1}{\partial x} = 0, \end{array} \right. \quad (33) $$

which, after introducing the new variables $u = \phi_1 + \phi_2$ and $v = (\phi_2 - \phi_1)/\sqrt{\alpha}$, can be rewritten as

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial u}{\partial t} + \frac{\partial v}{\partial x} = 0 \\[1em] \displaystyle v = \alpha \frac{\partial u}{\partial x}. \end{array} \right. \quad (34) $$

Once again, we recover the continuity equation and Fick's law, albeit with the original scalar field expressed as a particular linear combination of density and flux components.

More generally, for higher-dimensional diffusion equations we find that writing

$$ \begin{array}{l} \displaystyle \left(\frac{\partial}{\partial t} - \alpha \left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\right)\right) \phi(x, y, t) = 0 \quad (2D) \\[1em] \displaystyle \left(\eta_0 \frac{\partial}{\partial t} + i\sqrt{\alpha} \left(\eta_1 \frac{\partial}{\partial x} + \eta_2 \frac{\partial}{\partial y}\right) + \eta_c\right)^2 \vec{\phi}(x, y, t) = \vec{0}, \end{array} \quad (35) $$

in two dimensions, and

$$ \begin{array}{l} \displaystyle \left(\frac{\partial}{\partial t} - \alpha \left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}\right)\right) \phi(x, y, z, t) = 0 \quad (36) \\[1em] \displaystyle \left(\eta_0 \frac{\partial}{\partial t} + i\sqrt{\alpha} \left(\eta_1 \frac{\partial}{\partial x} + \eta_2 \frac{\partial}{\partial y} + \eta_3 \frac{\partial}{\partial z}\right) + \eta_c\right)^2 \vec{\phi}(x, y, z, t) = \vec{0}, \end{array} $$
in three dimensions, leads to the general set of conditions

$$ \eta_0^2 = 0, \quad \eta_i^2 = 1, \quad \eta_c^2 = 0, \quad \{\eta_0, \eta_i\} = 0, \quad \{\eta_0, \eta_c\} = 1, \quad \{\eta_i, \eta_c\} = 0, \quad (37) $$

where $i = 1, 2, \dots, d$. A minimal set of matrices satisfying these relations can be constructed as

$$ \eta_0 = \frac{e_{n+1} + ie_{n+2}}{2}, \quad \eta_i = e_i, \quad \eta_c = \frac{e_{n+1} - ie_{n+2}}{2}, \quad (38) $$

where $e_1, \dots, e_{n+2}$ denote matrix representations of the Clifford algebra $Cl_{n+2,0}(\mathbb{C})$ (or of higher-order Clifford algebras). We note that, due to the algebraic structure of these matrices, we cannot isolate the time-derivative operator. Consequently, we cannot write a closed time-evolution equation for the full state vector $\vec{\phi}$; instead, only a subset of its components is dynamically independent, while the remaining components are constrained algebraically without involving time derivatives.

Although this factorization generally yields systems with complex coefficients, we can choose matrix representations such that all resulting equations are real. In the case of the two-dimensional diffusion equation, this requires three real-valued matrices and one complex matrix, allowing the imaginary unit to be factored out explicitly. The standard set of four Dirac gamma matrices does not satisfy this requirement; however, by adopting alternative bases of $Cl_{5,0}(\mathbb{C})$, we can readily select suitable matrices. For example, taking $\gamma_i = \sigma_1 \otimes \sigma_i$ ($i = 1, 2, 3$) together with $\gamma_4 = \sigma_2 \otimes 1$ and $\gamma_5 = \sigma_3 \otimes 1$, and choosing the index set $(1, 3, 5, 2)$, we find, upon setting $\alpha = 1$ for simplicity,

$$ \left\{ \begin{aligned} \frac{\partial \phi_1}{\partial y} + \frac{\partial \phi_3}{\partial x} - \phi_4 &= 0 \\ \frac{\partial \phi_2}{\partial y} + \frac{\partial \phi_3}{\partial t} - \frac{\partial \phi_4}{\partial x} &= 0 \\ \frac{\partial \phi_1}{\partial x} - \phi_2 - \frac{\partial \phi_3}{\partial y} &= 0 \\ \frac{\partial \phi_1}{\partial t} - \frac{\partial \phi_2}{\partial x} - \frac{\partial \phi_4}{\partial y} &= 0. \end{aligned} \right. \quad (39) $$

This structure suggests that identifying $(\phi_1, \phi_2, \phi_3, \phi_4) \to (\phi_1, -j_1, \phi_2, -j_2)$ yields the following generalized continuity equations and generalized Fick-type laws for two-dimensional diffusive processes:

$$ \left\{ \begin{aligned} \frac{\partial \phi_1}{\partial t} + \frac{\partial j_1}{\partial x} + \frac{\partial j_2}{\partial y} &= 0 \\ \frac{\partial \phi_2}{\partial t} - \frac{\partial j_1}{\partial y} + \frac{\partial j_2}{\partial x} &= 0 \\ j_1 &= -\frac{\partial \phi_1}{\partial x} + \frac{\partial \phi_2}{\partial y} \\ j_2 &= -\frac{\partial \phi_1}{\partial y} - \frac{\partial \phi_2}{\partial x}. \end{aligned} \right. \quad (40) $$

Before turning to the next example, we emphasize that this type of factorization is not entirely new. Indeed, after an appropriate rescaling of variables, the factorization of the three-dimensional diffusion equation is structurally equivalent to the Lévy–Leblond
equation, which provides a Dirac-type linearization of the Schrödinger equation, or, equivalently, of the non-relativistic energy-momentum relation [16, 22]. Nevertheless, we stress that different choices of matrices, different identifications of the state-vector components, and different scalings may still lead to distinct physical interpretations.

### Example 3: Higher-Order and Fractional Roots

As a next example, we examine how the Dirac factorization method can be extended to higher-order roots and to non-integer, that is, fractional, derivatives. Up to this point, we have restricted ourselves to constructing square roots of operators. However, we can equally well consider higher roots, such as the cube root of an operator. To this end, we consider the simple pseudo-differential operator

$$ \sqrt[3]{1 + \alpha \frac{\partial^3}{\partial x^3}} \phi(x). \qquad (41) $$

There exist several ways to interpret such an operator. For instance, one may employ Laplace transform techniques to obtain a representation of the operator in the form

$$ \sqrt[3]{1 + \alpha \frac{\partial^3}{\partial x^3}} \phi(x) = \alpha^{-\frac{2}{3}} \int_0^\infty d\xi \, \xi_0 F_2\left(0; 1, \frac{4}{3}; -\frac{\xi^3}{27\alpha}\right) \left(1 + \alpha \frac{\partial^3}{\partial x^3}\right) \phi(x - \xi), \quad (42) $$

where $_0F_2$ denotes a generalized hypergeometric function. While this representation is formally well defined, it may be inconvenient in practice due to its intrinsic non-locality or its reliance on strong growth conditions imposed on the function $\phi$. It can therefore be advantageous to express the pseudo-differential operator without invoking an integral kernel. In close analogy to the standard Dirac factorization, this task reduces to expressing the cube root in the form [3]

$$ \sqrt[3]{a^3 + b^3} \mathbf{1} = a\tau_1 + b\tau_2, \qquad (43) $$

where $\tau_1$ and $\tau_2$ are generalized Pauli-Dirac matrices that satisfy the conditions

$$ \tau_i^2 = \mathbf{1}, \quad \tau_i^2 \tau_j + \tau_j \tau_i \tau_j + \tau_j \tau_i^2 = \mathbf{0}, \quad i, j \in \{1, 2\} \qquad (44) $$

The requirement that the cube of each matrix equals the identity implies that their eigenvalues must be powers of $\omega_3$, the cubic root of unity, which obeys the relations

$$ \omega_3 = e^{i\frac{2\pi}{3}}, \quad 1 + \omega_3 + \omega_3^2 = 0, \quad 1 \cdot \omega_3 \cdot \omega_3^2 = 1 \qquad (45) $$

Exploiting the properties of $\omega_3$, we can construct an ansatz for the first matrix as

$$ \tau_1 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \omega_3 & 0 \\ 0 & 0 & \omega_3^2 \end{pmatrix} \qquad (46) $$

To construct the second matrix, we consider the permutation of the ordered triple $(1 \ 2 \ 3) \to (2 \ 3 \ 1)$. The matrix indexing this permutation is known as a permutation matrix. Starting from the identity, the corresponding permutation matrix is given by

$$ \mathbf{P} = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} \qquad (47) $$
Since applying this permutation three times restores the original ordering, the permutation matrix satisfies, by construction,

$$ \mathbf{P}^3 = \mathbf{1}. \qquad (48) $$

One can readily verify that this matrix also fulfills the remaining algebraic conditions, thereby identifying $\mathbf{P}$ as the desired second matrix $\tau_2$. With this pair of matrices, we can express the pseudo-differential operator in the form [3]

$$ \sqrt[3]{1 + \alpha \frac{\partial^3}{\partial x^3}} \vec{\phi}(x) = \left( \tau_1 + \sqrt[3]{\alpha} \tau_2 \frac{\partial}{\partial x} \right) \vec{\phi}(x) = \begin{cases} \phi_1 + \sqrt[3]{\alpha} \frac{\partial \phi_2}{\partial x} \\ \omega_3 \phi_2 + \sqrt[3]{\alpha} \frac{\partial \phi_3}{\partial x} \\ \sqrt[3]{\alpha} \frac{\partial \phi_1}{\partial x} + \omega_3^2 \phi_3, \end{cases} \qquad (49) $$

where we have promoted the scalar function $\phi$ to a three-component vector $\vec{\phi}$. Again, this enlargement of the state space is the price to be paid to eliminate the non-local pseudo-differential operator. With this particular parametrization, all components of the vector-valued function are coupled. However, we can introduce a third matrix,

$$ \tau_3 = \frac{i}{\sqrt{3}} (\tau_2 \tau_1 - \tau_1 \tau_2) = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & \omega_3 \\ \omega_3^2 & 0 & 0 \end{pmatrix} \qquad (50) $$

which possesses the same algebraic properties as the other matrices. This allows us to construct a factorization that acts only on the individual components themselves:

$$ \sqrt[3]{1 + \alpha \frac{\partial^3}{\partial x^3}} \vec{\phi}(x) = \left( \tau_3 + \sqrt[3]{\alpha} \tau_2 \frac{\partial}{\partial x} \right) \vec{\phi}(x) = \begin{cases} \phi_2 + \sqrt[3]{\alpha} \frac{\partial \phi_2}{\partial x} \\ \omega_3^2 \phi_3 + \sqrt[3]{\alpha} \frac{\partial \phi_3}{\partial x} \\ \omega_3 \phi_1 + \sqrt[3]{\alpha} \frac{\partial \phi_1}{\partial x}. \end{cases} \qquad (51) $$

Employing this matrix-valued factorization of the pseudo-differential operator, the time-evolution equation

$$ \frac{\partial \phi}{\partial t} = \sqrt[3]{1 + \alpha \frac{\partial^3}{\partial x^3}} \phi(x) \qquad (52) $$

is transformed into a system of linear partial differential equations,

$$ \frac{\partial \vec{\phi}}{\partial t} = \left( \tau_i + \sqrt[3]{\alpha} \tau_j \frac{\partial}{\partial x} \right) \vec{\phi}(x) \qquad (53) $$

at the expense of introducing a vector-valued state function. The corresponding solution can then be written straightforwardly using the theory of time-evolution operators as

$$ \vec{\phi}(x, t) = \exp\left(t\left(\tau_i + \sqrt[3]{\alpha} \tau_j \frac{\partial}{\partial x}\right)\right) \vec{\phi}(x, 0), \qquad (54) $$

where the explicit form of the matrix exponential is provided in Appendix A for the appropriate choice $a_i = 1$ and $a_j = \sqrt[3]{\alpha} \frac{\partial}{\partial x}$.
We now turn to an example involving a fractional differential operator. Consider the pseudo-differential operator

$$ \sqrt{1 + \alpha \frac{\partial}{\partial x}} \phi(x). \qquad (55) $$

In this context, non-integer differential operators are commonly invoked to deal with the expression. However, this is not strictly necessary. A direct factorization using Pauli matrices yields

$$ \sqrt{1 + \alpha \frac{\partial}{\partial x}} \vec{\phi}(x) = \left( \sigma_i + \sqrt{\alpha} \sigma_j \hat{D}_x^{\frac{1}{2}} \right) \vec{\phi}(x), \qquad (56) $$

where $\hat{D}_x^{\frac{1}{2}}$ denotes the half-order derivative satisfying $\hat{D}_x^{\frac{1}{2}}(\hat{D}_x^{\frac{1}{2}}\phi(x)) = \frac{d}{dx}\phi(x)$ and which can be defined through various equivalent constructions [15, 19, 17]. Alternatively, we can preserve a linear, non-fractional derivative by seeking more general factorizations of the form

$$ \sqrt{1 + \alpha \frac{\partial}{\partial x}} \vec{\phi}(x) = \left( \varphi_1 + \alpha \varphi_2 \frac{\partial}{\partial x} \right) \vec{\phi}(x), \qquad (57) $$

which leads to the algebraic conditions

$$ \varphi_1^2 = 1, \quad \varphi_2^2 = 0, \quad \{\varphi_1, \varphi_2\} = 1 \qquad (58) $$

One set of matrices satisfying these relations is given, for example, by

$$ \varphi_1 = \sigma_i, \quad \varphi_2 = \frac{\sigma_i \pm i\sigma_j}{2}. \qquad (59) $$

In particular, choosing $(i, j) = (1, 2)$ in the first approach yields

$$ \sqrt{1 + \alpha \frac{\partial}{\partial x}} \vec{\phi}(x) = \begin{cases} \phi_2 - i\sqrt{\alpha} \hat{D}_x^{\frac{1}{2}} \phi_2 \\ \phi_1 + i\sqrt{\alpha} \hat{D}_x^{\frac{1}{2}} \phi_1, \end{cases} \qquad (60) $$

whereas setting $(i, j) = (3, 2)$ in the second approach gives

$$ \sqrt{1 + \alpha \frac{\partial}{\partial x}} \vec{\phi}(x) = \begin{cases} \phi_1 + \frac{\alpha}{2} \frac{\partial \phi_1}{\partial x} - \frac{\alpha}{2} \frac{\partial \phi_2}{\partial x} \\ \frac{\alpha}{2} \frac{\partial \phi_1}{\partial x} - \phi_2 - \frac{\alpha}{2} \frac{\partial \phi_2}{\partial x}. \end{cases} \qquad (61) $$

To illustrate how such an approach, which expresses operator roots through more elaborate factorizations rather than by naively taking roots of individual components, can uncover hidden physical structure, we consider the fractional pseudo-differential operator $\sqrt{1 + \alpha \frac{\partial^3}{\partial t^3}} \phi(t)$. A conventional factorization would lead to a $\frac{3}{2}$-order fractional time derivative via

$$ \sqrt{1 + \alpha \frac{\partial^3}{\partial t^3}} \vec{\phi}(t) = \left( \sigma_i + \sqrt{\alpha} \sigma_j \hat{D}_t^{\frac{3}{2}} \right) \vec{\phi}(t), \qquad (62) $$

However, we instead consider the ansatz

$$ \sqrt{1 + \alpha \frac{\partial^3}{\partial t^3}} \vec{\phi}(t) = \left( \boldsymbol{\mu}_1 + \boldsymbol{\mu}_2 \frac{\partial}{\partial t} + \boldsymbol{\mu}_3 \frac{\partial^2}{\partial t^2} \right) \vec{\phi}(t), \qquad (63) $$
which we interpret as a “damped harmonic oscillator”-type differential operator. Squaring this expression and expanding the resulting terms yields the set of algebraic conditions

$$ \mu_1^2 = 1, \quad \mu_2^2 + \{\mu_1, \mu_3\} = 0, \quad \mu_3^2 = 0, \quad \{\mu_1, \mu_2\} = 0, \quad \{\mu_2, \mu_3\} = \alpha \mathbf{1} \quad (64) $$

One particular choice satisfying these relations is

$$ \mu_1 = \sigma_i, \quad \mu_2 = -i \frac{\alpha}{2} \sigma_j + \frac{\sqrt{\alpha^2 - 8}}{2} \sigma_k, \quad \mu_3 = \sigma_i + i \sigma_j \quad (65) $$

Upon further choosing $(i, j, k) = (1, 2, 3)$, we obtain

$$ \sqrt{1 + \alpha \frac{\partial^3}{\partial t^3}} \vec{\phi}(t) = \begin{cases} \frac{\sqrt{\alpha^2 - 8}}{2} \frac{\partial \phi_1}{\partial t} + \phi_2 - \frac{\alpha}{2} \frac{\partial \phi_2}{\partial t} + 2 \frac{\partial^2 \phi_2}{\partial t^2} & \phi_1 + \frac{\alpha}{2} \frac{\partial \phi_1}{\partial t} - \frac{\sqrt{\alpha^2 - 8}}{2} \frac{\partial \phi_2}{\partial t} = 0 \\ \frac{\sqrt{\alpha^2 - 8}}{2} \frac{\partial \phi_1}{\partial t} + \phi_2 - \frac{\alpha}{2} \frac{\partial \phi_2}{\partial t} + 2 \frac{\partial^2 \phi_2}{\partial t^2} & \phi_1 + \frac{\alpha}{2} \frac{\partial \phi_1}{\partial t} - \frac{\sqrt{\alpha^2 - 8}}{2} \frac{\partial \phi_2}{\partial t} = 0 \end{cases} \quad (66) $$

This result demonstrates how the time-fractional pseudo-differential operator can be reinterpreted as describing two coupled components of a harmonic oscillator with effective damping or excitation, depending on the sign of $\alpha$, together with an equation containing a source term. The appearance of effective damping or excitation is consistent with analyses performed in the context of fractional harmonic oscillators of the form $(\hat{D}_t^\nu + k)\phi(t) = 0$, where one finds an effective friction term for $0 < \nu < 2$ and an excitation term for $2 < \nu$ [15], highlighting how these approaches, while distinct, lead to similar conclusions.

## Example 4: Telegrapher Equation

A particularly illustrative example is provided by the telegrapher equation. This equation describes a wide range of relativistic transport phenomena, and its solutions are characterized by a finite propagation speed. Phenomenologically, it is equivalent either to a wave equation supplemented by a dissipative term [1], modeled through a first-order time derivative, or to a diffusion equation in which the flux is governed by $\frac{\partial j}{\partial t} + j = \frac{\partial \phi}{\partial x}$ instead of the usual constitutive relation $j = \frac{\partial \phi}{\partial x}$ [4]. From a mathematical perspective, the telegrapher equation constitutes a hyperbolic differential equation and can be written as

$$ \frac{\partial^2 \phi}{\partial t^2} + 2 \frac{\partial \phi}{\partial t} - \frac{\partial^2 \phi}{\partial x^2} = 0. \quad (67) $$

Any choice of constant coefficients can be absorbed by an appropriate rescaling of space and time; consequently, we restrict our attention to this representative form. Taking the root of this equation provides a natural opportunity to explore different strategies for treating higher-order derivatives, each of which carry distinct physical implications.

As a first possibility, we may directly apply the standard Dirac factorization and, in the minimal setting, employ Pauli matrices to obtain [11]

$$ (I) \qquad \left( \sigma_i \frac{\partial}{\partial t} + \sqrt{2} \sigma_j \hat{D}_t^{\frac{1}{2}} + i \sigma_k \frac{\partial}{\partial x} \right) \vec{\phi}(x, t) = \vec{0}, \quad (68) $$

where a $\frac{1}{2}$-order fractional time derivative is introduced. Owing to the intrinsically non-local nature of fractional-order time derivatives, this linearization retains memory of the
past states of $\phi$; that is, it explicitly exhibits temporal memory effects. Alternatively, we may rewrite the original telegrapher equation in the equivalent form

$$ \left[ \left( \frac{\partial}{\partial t} + 1 \right)^2 - \left( \frac{\partial^2}{\partial x^2} + 1 \right) \right] \phi(x, t) = 0. \qquad (69) $$

Applying the usual minimal Dirac factorization to this expression yields an equation involving a pseudo-differential operator in space,

$$ (II) \qquad \left( \boldsymbol{\sigma}_i \frac{\partial}{\partial t} + \boldsymbol{\sigma}_i + i \boldsymbol{\sigma}_j \sqrt{1 + \frac{\partial^2}{\partial x^2}} \right) \vec{\phi}(x, t) = \vec{0}. \qquad (70) $$

In this case, the pseudo-differential operator $\sqrt{1 + \frac{\partial^2}{\partial x^2}} \phi$, which is typically defined through a convolution with a non-local integral kernel, induces spatial memory effects. This behavior stands in contrast to factorization (I), where the memory effects are purely temporal. Put differently, while the first factorization retains information about the past history of $\phi(x, t)$, the second factorization (II) yields a time evolution that depends on spatial information from all other points of the function at a given instant.

These are not the only possible factorizations. We can also identify a root of the telegrapher equation that remains local in both space and time and depends solely on first-order derivatives, namely

$$ (III) \qquad \left( \boldsymbol{\sigma}_i \frac{\partial}{\partial t} + i \boldsymbol{\sigma}_j \frac{\partial}{\partial x} + \boldsymbol{\sigma}_i + i \boldsymbol{\sigma}_k \right) \vec{\phi}(x, t) = \vec{0}. \qquad (71) $$

A solution to equation (III) follows straightforwardly from the theory of evolution operators. Using the identity $\boldsymbol{\sigma}_i^{-1} = \boldsymbol{\sigma}_i$, we rewrite the equation as

$$ \frac{\partial \vec{\phi}}{\partial t} = - \left( i \boldsymbol{\sigma}_i \boldsymbol{\sigma}_j \frac{\partial}{\partial x} + i \boldsymbol{\sigma}_i \boldsymbol{\sigma}_k + 1 \right) \vec{\phi}(x, t) \qquad (72) $$

and consequently obtain, for a given initial state $\vec{\phi}(x, 0)$,

$$ \vec{\phi}(x, t) = \exp \left( -t \left( i \boldsymbol{\sigma}_i \boldsymbol{\sigma}_j \frac{\partial}{\partial x} + i \boldsymbol{\sigma}_i \boldsymbol{\sigma}_k + 1 \right) \right) \vec{\phi}(x, 0). \qquad (73) $$

This expression can be treated using standard techniques, such as the Zassenhaus formula, or by exploiting the exponential of a Pauli vector,

$$ e^{it(a_1 \boldsymbol{\sigma}_1 + a_2 \boldsymbol{\sigma}_2 + a_3 \boldsymbol{\sigma}_3)} = \cos(\sqrt{a_1^2 + a_2^2 + a_3^2} t) \mathbf{1} + i \frac{\sin(\sqrt{a_1^2 + a_2^2 + a_3^2} t)}{\sqrt{a_1^2 + a_2^2 + a_3^2}} (a_1 \boldsymbol{\sigma}_1 + a_2 \boldsymbol{\sigma}_2 + a_3 \boldsymbol{\sigma}_3) \qquad (74) $$

with the specific choices $a_\ell = -i\epsilon_{ij\ell} \frac{\partial}{\partial x}$ and $a_m = -i\epsilon_{ikm}$, where $\epsilon$ denotes the Levi-Civita symbol. This procedure yields

$$ \vec{\phi}(x, t) = e^{-t} \left[ \cosh \left( t \sqrt{1 + \frac{\partial^2}{\partial x^2}} \right) \mathbf{1} - \frac{\sinh \left( t \sqrt{1 + \frac{\partial^2}{\partial x^2}} \right)}{\sqrt{1 + \frac{\partial^2}{\partial x^2}}} \left( \epsilon_{ij\ell} \boldsymbol{\sigma}_\ell \frac{\partial}{\partial x} + \epsilon_{ikm} \boldsymbol{\sigma}_m \right) \right] \vec{\phi}(x, 0), \qquad (75) $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Two solutions to the unconventional root (IV) of the telegrapher equation for a Gaussian initial condition (left) and a box initial condition (right).</figcaption></figure>

where we have assumed $i \neq j \neq k$.

Furthermore, one may consider even more unconventional factorizations, such as

$$ (IV) \qquad \left( \sqrt{2}\sigma_i \hat{D}_t^{\frac{1}{2}} + \sigma_j \sqrt{\frac{\partial^2}{\partial t^2} - \frac{\partial^2}{\partial x^2}} \right) \vec{\phi}(x,t) = \vec{0}, \qquad (76) $$

which combines a fractional time derivative with a pseudo-differential operator given by the square root of the wave operator. This root is neither local in time nor in space and is best understood using a combined Fourier–Laplace transform,

$$ (\sigma_i \sqrt{2s} + \sigma_j \sqrt{s^2 + k^2}) \vec{\Phi}(s, k) = \vec{0}, \qquad (77) $$

where

$$ \vec{\Phi}(s, k) = \int_{0}^{\infty} dt \int_{-\infty}^{\infty} dx \, e^{-st-ikx} \vec{\phi}(x, t). \qquad (78) $$

For the specific choice $(i, j) = (1, 2)$, the equations decouple, and we illustrate the corresponding solutions for a Gaussian and box initial condition in Fig. 3.

### Example 5: Dynamic Euler–Bernoulli Equation

The dynamics of thin beams undergoing small deflections in the absence of distributed loads are governed by the fourth-order Euler–Bernoulli equation [12]

$$ \left( \frac{\partial^2}{\partial t^2} + \alpha \frac{\partial^4}{\partial x^4} \right) w(x, t) = 0, \qquad (79) $$

where $w$ denotes the transverse deflection of the beam perpendicular to the $x$-direction. Depending on the imposed boundary conditions, the solutions exhibit vibratory motion
that partially mimics wave propagation in an inhomogeneous medium. While we may factor this equation by directly taking its root, which leads to two “diffusion”-type (or, depending on the choice of indices, also “Schrödinger”-type) equations [6],

$$ \left( \sigma_i \frac{\partial}{\partial t} + \sqrt{\alpha} \sigma_j \frac{\partial^2}{\partial x^2} \right)^2 \vec{w}(x, t) = \vec{0}, \qquad (80) $$

the wave-like nature of the original equation suggests looking for a more explicitly “wave”-type factorizations of the form

$$ \left( \beta_1 \frac{\partial^2}{\partial t^2} + \sqrt{\alpha} \beta_2 \frac{\partial^2}{\partial x^2} + \beta_3 \right)^2 \vec{w}(x, t) = \vec{0}, \qquad (81) $$

in which both temporal and spatial derivatives appear at second order. The associated matrices are algebraically analogous to those arising in the factorization of the diffusion equation and can therefore be chosen, in the simplest case, as

$$ \beta_1 = \frac{\sigma_i \pm i\sigma_j}{2}, \quad \beta_2 = \sigma_k, \quad \beta_3 = \frac{\sigma_i \mp i\sigma_j}{2}. \qquad (82) $$

A particular choice $(i, j, k) = (1, 2, 3)$ yields the coupled system

$$ \left\{ \begin{array}{l} \frac{\partial^2 w_1}{\partial t^2} - \sqrt{\alpha} \frac{\partial^2 w_2}{\partial x^2} = 0 \\ w_2 + \sqrt{\alpha} \frac{\partial^2 w_1}{\partial x^2} = 0, \end{array} \right. \qquad (83) $$

which can be further rearranged into

$$ \left\{ \begin{array}{l} \frac{\partial^2 w_2}{\partial t^2} - \sqrt{\alpha} \frac{\partial^2 w_2}{\partial x^2} = w_1 + \sqrt{\alpha} \frac{\partial^2 w_1}{\partial x^2} \\ \frac{\partial^2 w_1}{\partial t^2} + \alpha \frac{\partial^4 w_1}{\partial x^4} = 0. \end{array} \right. \qquad (84) $$

In this form, the wave character of $w_2$ becomes explicit, while $w_1$ satisfies an inhomogeneous wave equation with a source term determined by $w_1$ itself through the dynamic Euler-Bernoulli equation. Since the physical behavior of the system is already fully encoded in $w_1$, this particular rearrangement does not provide additional physical insight.

A more compelling choice for the matrices $\beta_1$, $\beta_2$, and $\beta_3$ is given by

$$ \beta_1 = \begin{pmatrix} 1 & 1 \\ -1 & -1 \end{pmatrix}, \quad \beta_2 = \begin{pmatrix} -1 & -2 \\ 0 & 1 \end{pmatrix}, \quad \beta_3 = \begin{pmatrix} 0 & -1 \\ 0 & 0 \end{pmatrix}. \qquad (85) $$

With this choice, we obtain the system

$$ \left\{ \begin{array}{l} \frac{\partial^2 w_1}{\partial t^2} - \sqrt{\alpha} \frac{\partial^2 w_1}{\partial x^2} = w_2 - \frac{\partial^2 w_2}{\partial t^2} + 2\sqrt{\alpha} \frac{\partial^2 w_2}{\partial x^2} \\ \frac{\partial^2 w_1}{\partial t^2} = -\frac{\partial^2 w_2}{\partial t^2} + \sqrt{\alpha} \frac{\partial^2 w_2}{\partial x^2}, \end{array} \right. \qquad (86) $$
which can be rearranged and, upon identifying $w_1 \to w$ as the beam deflection and $w_2 \to \sqrt{\alpha}s$ as a source or inhomogeneity, yields

$$ \left\{ \begin{array}{l} \frac{2}{\sqrt{\alpha}} \frac{\partial^2 w}{\partial t^2} - \frac{\partial^2 w}{\partial x^2} = s - 2 \frac{\partial^2 s}{\partial t^2} + 3 \sqrt{\alpha} \frac{\partial^2 s}{\partial x^2} \\ \\ \sqrt{\alpha} \frac{\partial^2 s}{\partial x^2} + s = - \frac{\partial^2 w}{\partial x^2}. \end{array} \right. \qquad (87) $$

This formulation indicates that, under the present factorization, the deflection $w$ is described by waves with propagation velocity $c = \sqrt{\frac{\sqrt{\alpha}}{2}}$, sourced by a field $s$ that is determined through an inhomogeneous Helmholtz equation with effective wave number $k = 1/\sqrt{\alpha}$ and driving term proportional to $w$.

Alternatively, we may factor the Euler–Bernoulli equation into a continuity equation and a generalized “Euler–Bernoulli–Fick” law by considering

$$ \left( \beta_1 \frac{\partial}{\partial t} + \alpha \beta_2 \frac{\partial^3}{\partial x^3} + \beta_3 \frac{\partial}{\partial x} \right)^2 \vec{\phi}, \qquad (88) $$

with matrix coefficients chosen as

$$ \beta_1 = \sigma_i, \quad \beta_2 = \frac{\sigma_j + i\sigma_k}{2}, \quad \beta_3 = \frac{\sigma_j - i\sigma_k}{2}. \qquad (89) $$

For the specific choice $(i, j, k) = (3, 1, 2)$, we obtain

$$ \left\{ \begin{array}{l} \frac{\partial w}{\partial t} + \frac{\partial j_w}{\partial x} = 0 \\ \\ \frac{\partial j_w}{\partial t} = \alpha \frac{\partial^3 w}{\partial x^3}, \end{array} \right. \qquad (90) $$

where $j_w$ may be interpreted as a “deflection flux.”

Finally, we can attempt a complete linearization of the dynamic Euler–Bernoulli equation by seeking matrices such that

$$ \left( \frac{\partial^2}{\partial t^2} + \frac{\partial^4}{\partial x^4} \right) \vec{w}(x, t) = \left( \sigma_i \frac{\partial}{\partial t} + \sigma_j \frac{\partial^2}{\partial x^2} \right)^2 \vec{w}(x, t), \qquad (91) $$

and subsequently expressing

$$ \left( \sigma_i \frac{\partial}{\partial t} + \sigma_j \frac{\partial^2}{\partial x^2} \right) \vec{w}(x, t) = \left( \delta_1 \frac{\partial}{\partial t} + \delta_2 \frac{\partial}{\partial x} + \delta_3 \right)^2 \vec{w}(x, t). \qquad (92) $$

This construction leads to the algebraic conditions

$$ \delta_1^2 = 0, \quad \delta_2^2 = \sigma_j, \quad \delta_3^2 = 0, \quad \{\delta_1, \delta_3\} = \sigma_i, \quad \{\delta_1, \delta_2\} = 0, \quad \{\delta_2, \delta_3\} = 0, \qquad (93) $$

which admit no solutions in the $2 \times 2$ case. However, by generalizing the Pauli matrices to $4 \times 4$ matrices through the construction $\mathbb{1}_2 \otimes \sigma_i$, we obtain the solutions

$$ \delta_1 = \mathbf{E} \otimes \sigma_i, \quad \delta_2 = \sigma_3 \otimes \sqrt{\sigma_j}, \quad \delta_3 = \mathbf{E}^\mathrm{T} \otimes \mathbb{1}_2, \qquad (94) $$
<figure><img src="image_5.png" /><figcaption>Figure 4: Plate immersed in a Newtonian fluid connected by a massless spring (figure adapted from [13])</figcaption></figure>

where we define

$$ \sqrt{\sigma_j} = \frac{1-i}{2}(\sigma_j + i\mathbb{1}_2), \quad \mathbf{E} = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}, \qquad (95) $$

with the useful properties

$$ \{\mathbf{E}, \mathbf{E}^T\} = \mathbb{1}_2, \quad \{\sigma_3, \mathbf{E}\} = \{\sigma_3, \mathbf{E}^T\} = \mathbb{0}_2. \qquad (96) $$

In this way, we can factor the original dynamic Euler-Bernoulli equation into the transport-type form

$$ \left(\frac{\partial^2}{\partial t^2} + \alpha \frac{\partial^4}{\partial x^4}\right) \ddot{w}(x,t) = \left(\delta_1 \frac{\partial}{\partial t} + \sqrt[4]{\alpha} \delta_2 \frac{\partial}{\partial x} + \delta_3\right)^4 \dot{w}(x,t). \qquad (97) $$

Unfortunately, none of the index choices yield a system of purely real equations, complicating their physical interpretation.

### Example 6: Free Bagley-Torvik Equation

The Bagley-Torvik equation describes the motion of a rigid plate of mass $m$ and area $A$ immersed in a Newtonian fluid of infinite extent with density $\rho$ and viscosity $\mu$ [21]. The plate is attached to a fixed point by a massless spring with stiffness $k$, see Fig. 4. Considering the displacement $h(t)$ relative to a fixed reference position $h_0$ and neglecting external forces, we find that the motion is governed by the fractional-order differential equation [21]

$$ \left(m \frac{d^2}{dt^2} + 2A \sqrt{\rho \mu} \hat{D}_t^{\frac{3}{2}} + k \right) h(t) = 0, \qquad (98) $$

where the fractional derivative is understood in the Caputo sense [8]. For notational convenience, we define $\alpha = \frac{2A\sqrt{\rho\mu}}{m}$ and $\beta = \frac{k}{m}$, which allows us to rewrite the equation in
<figure><img src="image_6.png" /><figcaption>Figure 5: Solutions to the higher-order representation of the Bagley–Torvik equation (102) for vanishing boundary conditions and different choices of parameters. In both cases $m = 1$ and $\rho\mu = 1$.</figcaption></figure>

the compact form

$$ \left( \frac{d^2}{dt^2} + \alpha \hat{D}_t^{\frac{3}{2}} + \beta \right) h(t) = 0. \qquad (99) $$

As discussed previously, taking roots of differential operators naturally gives rise to fractional operators. Viewed conversely, we may interpret a fractional-order differential equation as representing a suitable root of a higher-order, non-fractional differential equation. One possible realization of this idea is to assume that the original equation corresponds to a single component (specifically $h_1$) of a vector-valued equation of the form

$$ \left( \sigma_1 \alpha \hat{D}_t^{\frac{3}{2}} + i \sigma_2 \left( \frac{d^2}{dt^2} + \beta \right) \right) \vec{h}(t) = \vec{0}. \qquad (100) $$

Squaring the operator recovers the associated “overlying” fourth-order differential equation,

$$ \left( \frac{d^4}{dt^4} - \alpha^2 \frac{d^3}{dt^3} + 2\beta \frac{d^2}{dt^2} + \beta^2 \right) \vec{h}(t) = \vec{0}. \qquad (101) $$

Consequently, all solutions of the original Bagley–Torvik equation are also solutions of this higher-order equation. We may rearrange the latter into a form reminiscent of Newton’s second law,

$$ m \frac{d^2 h}{dt^2} = -\frac{k}{2} h(t) + \frac{2A^2 \rho \mu}{k} \frac{d^3 h}{dt^3} - \frac{m^2}{2k} \frac{d^4 h}{dt^4}. \qquad (102) $$

In this representation, the system appears as a harmonic oscillator with effective spring constant $k/2$, in addition to a dissipative friction force that depends on higher-order time derivatives, $F_{\text{fric}} = -\frac{2A^2 \rho \mu}{k} \frac{d^3 h}{dt^3} + \frac{m^2}{2k} \frac{d^4 h}{dt^4}$. We illustrate the corresponding solutions and the associated friction force in Fig. 5 for several choices of parameters. It should be noted that this equation also admits additional solutions that do not satisfy the Bagley–Torvik equation and can only be distinguished through the specification of appropriate boundary conditions.
<figure><img src="image_7.png" /><figcaption>Figure 6: Solutions to the root system of the higher-order representation of the Bagley–Torvik equation (107) for vanishing boundary conditions and different choices of parameters. In both cases $m = 1$ and $\rho\mu = 1$.</figcaption></figure>

We may also seek a formulation of this higher-order differential equation in terms of non-fractional orders. To this end, we repeat the factorization procedure and assume a more tractable, explicit “Newtonian” form,

$$ \left( \kappa_1 \frac{d^2}{dt^2} + \kappa_2 \frac{d}{dt} + \kappa_3 \right)^2 \vec{h}(t) = \vec{0}. \qquad (103) $$

From this ansatz, we obtain the algebraic conditions

$$ \kappa_1^2 = 1, \quad \kappa_2^2 + \{\kappa_1, \kappa_3\} = 2\beta 1, \quad \kappa_3^2 = \beta^2 1, \quad \{\kappa_1, \kappa_2\} = -\alpha^2 1, \quad \{\kappa_2, \kappa_3\} = 0. \quad (104) $$

One particular solution of these constraints is

$$ \kappa_1 = \sigma_i, \quad \kappa_2 = -\frac{\alpha^2}{2}\sigma_i + i\frac{\sqrt{\alpha^4 - 8\beta}}{2}\sigma_j, \quad \kappa_3 = \beta\sigma_k. \qquad (105) $$

Rearranging the resulting equations yields a form closely resembling Newton’s second law,

$$ m \frac{d^2 \vec{h}}{dt^2} = -m \kappa_1 \kappa_2 \frac{d \vec{h}}{dt} - m \kappa_1 \kappa_3 \vec{h}. \qquad (106) $$

Here, the right-hand side represents the total force acting on the plate described by the state vector $\vec{h}$. The first term corresponds to a Stokes-type friction force with a highly non-trivial effective friction coefficient, while the second term is simply Hooke’s law, recalling that $h$ is measured relative to equilibrium. For the choice $(i, j, k) = (1, 2, 3)$, we obtain the comparatively simple system

$$ \left\{ \begin{aligned} m \frac{d^2 h_1}{dt^2} &= \frac{2}{m} \left( A^2 \rho \mu + \sqrt{(A^2 \rho \mu)^2 - \frac{k}{2m^3}} \right) \frac{dh_1}{dt} + kh_2, \\ m \frac{d^2 h_2}{dt^2} &= \frac{2}{m} \left( A^2 \rho \mu - \sqrt{(A^2 \rho \mu)^2 - \frac{k}{2m^3}} \right) \frac{dh_2}{dt} - kh_1, \end{aligned} \right. \qquad (107) $$

which admits decaying solutions and, notably, involves no fractional derivatives. We present a representative set of such solutions in Fig. 6.
## Final Remarks

In this short note, we have summarized a variety of different factorizations, in the spirit of Dirac, applied to common differential equations arising in mathematical physics. We have demonstrated that employing more elaborate forms of root operators can reveal novel insights into the underlying structure of these equations and have provided several concrete examples illustrating how the resulting phenomenology may be interpreted and understood. We showed that different ways of factoring partial differential equations, together with the multitude of possible combinations and arrangements of operators, can lead to fundamentally different manifestations of hidden structures. This behavior was exemplified by the telegrapher equation, for which one factorization exhibits time-memory effects, another retains spatial memory, while still other factorizations are fully linear in the derivative operators and thus give rise to entirely local linear equations.

Although the appearance of fractional derivative operators or pseudo-differential operators may at times seem unavoidable within the framework of standard Dirac factorizations, we have shown that this outcome can often be circumvented by adopting more intricate forms of the root operators. In this manner, we demonstrated how a well-known fractional differential equation – the Bagley–Torvik equation describing the motion of plates immersed in Newtonian fluids – can be reinterpreted as a non-fractional ordinary differential equation and subsequently factorized further to yield a system of Newtonian second-law equations with non-trivial friction coefficients.

This work has touched upon several distinct and, in part, loosely connected topics, many of which we have not explored in full detail or provided exhaustive interpretations for. Nevertheless, we believe that the general approach outlined in this note may prove fruitful for uncovering the underlying physics and mathematical structures of a wide range of natural processes and therefore merits further and more systematic investigation.

# Appendix

## A Matrix Exponential of Cube-Root Matrices $\tau$

To compute matrix exponentials involving the cube-root matrices $\tau$, specifically expressions of the form $e^{t(a_1\tau_1 + a_2\tau_2 + \dots)}$, we begin by writing the series definition of the exponential:

$$e^{t(\vec{a} \cdot \tau)} = \sum_{n=0}^{\infty} \frac{t^n}{n!} (\vec{a} \cdot \tau)^n, \qquad (A.1)$$

where we have introduced the cube-root matrix vector $(\vec{a} \cdot \tau) = a_1\tau_1 + a_2\tau_2 + \dots$. By definition, the relation $(\vec{a} \cdot \tau)^3n = (a_1^3 + a_2^3 + \dots)^n = \|\vec{a}\|_3^{3n}$ holds, which allows us to decompose the sum into three distinct contributions:

$$\begin{aligned} \sum_{n=0}^{\infty} \frac{t^n}{n!} (\vec{a} \cdot \tau)^n &= \sum_{n=0}^{\infty} \frac{t^{3n}}{(3n)!} (\vec{a} \cdot \tau)^{3n} + \sum_{n=0}^{\infty} \frac{t^{3n+1}}{(3n+1)!} (\vec{a} \cdot \tau)^{3n+1} + \sum_{n=0}^{\infty} \frac{t^{3n+2}}{(3n+2)!} (\vec{a} \cdot \tau)^{3n+2} \\ &= \sum_{n=0}^{\infty} \frac{(t\|\vec{a}\|_3)^{3n}}{(3n)!} \mathbb{1} + \sum_{n=0}^{\infty} \frac{(t\|\vec{a}\|_3)^{3n+1}}{(3n+1)!} \frac{(\vec{a} \cdot \tau)}{\|\vec{a}\|_3} + \sum_{n=0}^{\infty} \frac{(t\|\vec{a}\|_3)^{3n+2}}{(3n+2)!} \frac{(\vec{a} \cdot \tau)^2}{\|\vec{a}\|_3^2}. \end{aligned}$$
These sums correspond to the well-known trisections of the exponential series, given explicitly by

$$T_1(z) = \sum_{n=0}^{\infty} \frac{z^{3n}}{(3n)!} = \frac{1}{3} \left( e^z + 2e^{-\frac{z}{2}} \cos \left( \frac{\sqrt{3}z}{2} \right) \right),$$

$$T_2(z) = \sum_{n=0}^{\infty} \frac{z^{3n+1}}{(3n+1)!} = \frac{1}{3} \left( e^z - 2e^{-\frac{z}{2}} \cos \left( \frac{\sqrt{3}z}{2} + \frac{\pi}{3} \right) \right),$$

$$T_3(z) = \sum_{n=0}^{\infty} \frac{z^{3n+2}}{(3n+2)!} = \frac{1}{3} \left( e^z - 2e^{-\frac{z}{2}} \cos \left( \frac{\sqrt{3}z}{2} - \frac{\pi}{3} \right) \right).$$

Consequently, we obtain the compact expression

$$e^{t(\vec{a} \cdot \vec{\tau})} = T_1(t\|\vec{a}\|_3)\mathbb{1} + T_2(t\|\vec{a}\|_3) \frac{(\vec{a} \cdot \vec{\tau})}{\|\vec{a}\|_3} + T_3(t\|\vec{a}\|_3) \frac{(\vec{a} \cdot \vec{\tau})^2}{\|\vec{a}\|_3^2}. \quad (A.2)$$

## B Finding Coefficients in a Clifford Expansion as a Gram Factorization Problem

We consider a set of elements $\mathbf{A}_i$ in a Clifford algebra expressed as a linear combination of generators $\mathbf{e}_j$:

$$\mathbf{A}_i = \sum_{j=1}^{n} c_j^{(i)} \mathbf{e}_j, \quad i = 1, \dots, m, \qquad (B.3)$$

where the generators satisfy the anticommutation relations $\{\mathbf{e}_j, \mathbf{e}_k\} = 2\delta_{jk}$. Given the constraints on the squares $\mathbf{A}_i^2 = b_i$ and the anticommutators $\{\mathbf{A}_i, \mathbf{A}_k\} = b_{ik}$, we determine the coefficients $c_j^{(i)}$ by first expanding the square of $\mathbf{A}_i$:

$$\mathbf{A}_i^2 = \left( \sum_{j=1}^{n} c_j^{(i)} \mathbf{e}_j \right)^2 = \sum_{j,k} c_j^{(i)} c_k^{(i)} \mathbf{e}_j \mathbf{e}_k = \sum_{j=1}^{n} (c_j^{(i)})^2 + \sum_{j \neq k} c_j^{(i)} c_k^{(i)} \mathbf{e}_j \mathbf{e}_k. \quad (B.4)$$

The Clifford relations imply that the cross terms $\mathbf{e}_j \mathbf{e}_k$ for $j \neq k$ anticommute and cancel under the symmetric summation, yielding the condition $\sum_{j=1}^{n} (c_j^{(i)})^2 = b_i$. Similarly, we expand the anticommutators:

$$\{\mathbf{A}_i, \mathbf{A}_k\} = \sum_{j,\ell} c_j^{(i)} c_\ell^{(k)} \{\mathbf{e}_j, \mathbf{e}_\ell\} = 2 \sum_{j=1}^{n} c_j^{(i)} c_j^{(k)} = b_{ik}, \quad (B.5)$$

which simplifies to the requirement $\sum_{j=1}^{n} c_j^{(i)} c_j^{(k)} = \frac{b_{ik}}{2}$.

We reformulate these conditions by defining an $n \times m$ matrix $\mathbf{C}$ with entries $C_{ij} = c_j^{(i)}$:

$$\mathbf{C} = \begin{pmatrix} c_1^{(1)} & c_2^{(1)} & \cdots \\ c_1^{(2)} & c_2^{(2)} & \cdots \\ \vdots & \vdots & \ddots \end{pmatrix}. \qquad (B.6)$$
The scalar constraints correspond to the Gram matrix factorization $CC^T = B$, where $B$ is a symmetric matrix defined by $B_{ii} = b_i$ and $B_{ik} = \frac{b_{ik}}{2}$ for $i \neq k$. To solve for $C$, we utilize the eigendecomposition of the symmetric matrix $B = U\Lambda U^T$, where $U$ is an orthogonal matrix and $\Lambda$ is a diagonal matrix of eigenvalues. We then construct the coefficient matrix as:

$$C = U\Lambda^{\frac{1}{2}}Q, \qquad (B.7)$$

where $Q$ represents an arbitrary orthogonal matrix. Each row of the resulting matrix $C$ provides the coefficients $c_j^{(i)}$ necessary for the expansion.

## References

[1] Y.M. Ali and L.C. Zhang. "Relativistic heat conduction". In: *International Journal of Heat and Mass Transfer* 48.12 (2005). ISSN: 0017-9310. DOI: 10.1016/j.ijheatmasstransfer.2005.02.003. URL: http://dx.doi.org/10.1016/j.ijheatmasstransfer.2005.02.003.

[2] D. Babusci and G. Dattoli. *The Dirac factorization method and the harmonic oscillator*. 2012. DOI: 10.48550/ARXIV.1203.2526. URL: https://arxiv.org/abs/1203.2526.

[3] D. Babusci et al. "A note on the extension of the Dirac method". In: *Applied Mathematics and Computation* 218.4 (2011). ISSN: 0096-3003. DOI: 10.1016/j.amc.2011.06.002. URL: http://dx.doi.org/10.1016/j.amc.2011.06.002.

[4] Carlo Cattaneo. "Sur une forme de l'equation de la chaleur eliminant la paradoxe d'une propagation instantantee". In: *Compt. Rendu* 247 (1958), pp. 431-433.

[5] Jean le Rond d'Alembert. "Recherches sur la courbe que forme une corde tendue mise en vibration". In: (1747).

[6] Giuseppe Dattoli, Emanuele Di Palma, and Alessandro Curcio. "Dirac Factorization, Partial/Ordinary Differential Equations and Fractional Calculus". In: *Symmetry* 17.11 (2025). ISSN: 2073-8994. DOI: 10.3390/sym17111984. URL: http://dx.doi.org/10.3390/sym17111984.

[7] Giuseppe Dattoli and Amalia Torre. "Root Operators and "Evolution" Equations". In: *Mathematics* 3.3 (2015). ISSN: 2227-7390. DOI: 10.3390/math3030690. URL: http://dx.doi.org/10.3390/math3030690.

[8] Kai Diethelm. "General theory of Caputo-type fractional differential equations". In: *Fractional Differential Equations*. De Gruyter, Feb. 2019, pp. 1-20. ISBN: 9783110571660. DOI: 10.1515/9783110571660-001. URL: http://dx.doi.org/10.1515/9783110571660-001.

[9] Paul Adrien Maurice Dirac. "Relativistic quantum mechanics". In: *Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character* 136.829 (1932). ISSN: 2053-9150. DOI: 10.1098/rspa.1932.0094. URL: http://dx.doi.org/10.1098/rspa.1932.0094.
[10] Paul Adrien Maurice Dirac. "The quantum theory of the electron". In: *Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character* 117.778 (1928). ISSN: 2053-9150. DOI: 10.1098/rspa.1928.0023. URL: http://dx.doi.org/10.1098/rspa.1928.0023.

[11] Zine El Abiddine Fellah et al. "Factorization à la Dirac Applied to Some Equations of Classical Physics". In: *Mathematics* 9.8 (2021). ISSN: 2227-7390. DOI: 10.3390/math9080899. URL: http://dx.doi.org/10.3390/math9080899.

[12] James M. Gere and Stephen P. Timoshenko. *Mechanics of Materials*. Springer US, 1991. ISBN: 9781489931245. DOI: 10.1007/978-1-4899-3124-5. URL: http://dx.doi.org/10.1007/978-1-4899-3124-5.

[13] Juan Luis González-Santander and Alexander Apelblat. *Analytical closed-form solution of the Bagley-Torvik equation*. 2025. DOI: 10.48550/ARXIV.2501.05972. URL: https://arxiv.org/abs/2501.05972.

[14] Walter Greiner. *Relativistic Quantum Mechanics*. Springer Berlin Heidelberg, 1995. ISBN: 9783642880827. DOI: 10.1007/978-3-642-88082-7. URL: http://dx.doi.org/10.1007/978-3-642-88082-7.

[15] Richard Herrmann. *Fractional Calculus: An Introduction for Physicists*. WORLD SCIENTIFIC, 2013. ISBN: 9789814551083. DOI: 10.1142/8934. URL: http://dx.doi.org/10.1142/8934.

[16] Jean-Marc Lévy-Leblond. "Nonrelativistic particles and wave equations". In: *Communications in Mathematical Physics* 6.4 (1967). ISSN: 1432-0916. DOI: 10.1007/bf01646020. URL: http://dx.doi.org/10.1007/BF01646020.

[17] Kenneth S Miller and Bertram Ross. "An introduction to the fractional calculus and fractional differential equations". In: (No Title) (1993).

[18] W. Pauli. "Zur Quantenmechanik des magnetischen Elektrons". In: *Zeitschrift für Physik* 43.9-10 (1927). ISSN: 0044-3328. DOI: 10.1007/bf01397326. URL: http://dx.doi.org/10.1007/BF01397326.

[19] Igor Podlubny. *Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications*. Vol. 198. elsevier, 1998.

[20] Erwin Schroedinger. "The factorization of the hypergeometric equation". In: (1999). DOI: 10.48550/ARXIV.PHYSICS/9910003.

[21] P. J. Torvik and R. L. Bagley. "On the Appearance of the Fractional Derivative in the Behavior of Real Materials". In: *Journal of Applied Mechanics* 51.2 (1984). ISSN: 1528-9036. DOI: 10.1115/1.3167615. URL: http://dx.doi.org/10.1115/1.3167615.

[22] James M Wilkes. "The Pauli and Lévy-Leblond equations, and the spin current density". In: *European Journal of Physics* 41.3 (Mar. 2020), p. 035402. ISSN: 1361-6404. DOI: 10.1088/1361-6404/ab7495. URL: http://dx.doi.org/10.1088/1361-6404/ab7495.