Spectral Properties of Stochastic Processes Possessing Finite Propagation Velocity

Abstract

This article investigates the spectral structure of the evolution operators associated with the statistical description of stochastic processes possessing finite propagation velocity. Generalized Poisson–Kac processes and Lévy walks are explicitly considered as paradigmatic examples of regular and anomalous dynamics. A generic spectral feature of these processes is the lower boundedness of the real part of the eigenvalue spectrum that corresponds to an upper limit of the spectral dispersion curve, physically expressing the relaxation rate of a disturbance as a function of the wave vector. We also analyze Generalized Poisson–Kac processes possessing a continuum of stochastic states parametrized with respect to the velocity. In this case, there is a critical value for the wave vector, above which the point spectrum ceases to exist, and the relaxation dynamics becomes controlled by the essential part of the spectrum. This model can be extended to the quantum case, and in fact, it represents a simple and clear example of a sub-quantum dynamics with hidden variables.

1. Introduction

The investigation of micro- and nanoscale physics [1] as well as the extension of statistical physical concepts to new phenomelogies involving active matter and living beings [2,3] have stimulated the development of more refined descriptions of stochastic processes, accounting for background and thermal or quantum fluctuations, aimed at deriving their statistical features and long-term, in some cases anomalous [4], scaling properties [5,6]. Experimental results on the motion of bacteria, amoebae, insects, and other classes of living beings indicate that the classical paradigm of Wiener fluctuations (the mathematical Brownian motion) does not always apply to the erratic kinematics of these entities [7].

Turning our attention to a completely different branch of physics, namely field theory, the physical constraints on the texture of the space–time fabric dictate the requirement for the finite propagation velocity as an essential condition in relativistic consistent models of fluctuations [8,9]. What is remarkable in this context is that the Lorentz covariance applies to the opposite ends of the length scale spectrum, i.e., at very short (quantum fluctuations at the particle level) [10] and very large length scales (cosmological models) [11].

But even for processes in the everyday human experience, such as those related to mass, momentum, and heat transport that represent the natural realm of investigation of thermodynamic and transport theories, the quest to overcome the classical paradoxes resulting from the use of stochastic equations with infinite propagation speed—unavoidable within the assumption of constitutive equations of the Fickian type (i.e., where the flux of the thermodynamic entities is proportional and opposite to their concentration gradients, such as in the Fick, Fourier, and Newton constitutive equations for mass, heat, and momentum “diffusive” fluxes)—is not only dictated by aesthetical or epistemological arguments supporting the internal coherence of these macroscopic and statistical theories [12], but stems from the resolution of physical problems such as those of solute transport through polymeric materials or heat transport in nanodevices [13,14].

A common denominator in all these issues can be found in the understanding of the emergent properties characterizing stochastic models that possess a finite propagation velocity and of their distinctive features, either as regards the local regularity of their trajectories or their collective statistical behavior. In this article, we focus on the latter issue. Specifically, we study the spectral structure of the evolution operators that propagate the probability densities of these models. Henceforth, we refer to these operators as Statistical Evolution Operators, or SEO. This is because the constraint of bounded velocity enforces specific and characteristic spectral features, distinguishing them from their classical and widely used counterparts, namely the Wiener and Wiener-driven fluctuations. We thoroughly investigate the fluctuation spectra of Generalized Poisson-Kac processes (GPK) [15,16,17,18] and Lévy Walks (LW) [19,20,21,22], as these processes represent two of the main classes of stochastic models characterized by bounded propagation velocity.

GPK processes originate from the need to define a sufficiently wide class of stochastic models that are able to interpret the equations of extended thermodynamics for physical processes far from equilibrium at the mesoscopic level [12,23]. This represents an extension of the original model proposed by Kac [15] of a stochastic process driven by the parity of a Poisson counting process, which has been the subject of intense investigation in the past as a prototype of a non-Markovian stochastic motion driven by bounded, dichotomous, and colored noise [24,25,26,27,28]. We choose these two classes of processes because they are subjected to a common and complete statistical description in terms of the partial probability densities parametrized with respect to the internal variables of the process (for GPK, the velocity direction; for LWs, the transition age). This was shown especially for LWs by the recent work of Fedotov [29] (see also [30,31,32,33]), and by the unifying theory of Extended Poisson–Kac processes developed by Giona and collaborators [34]. Our main result is that the real part of the eigenvalues of the associated SEO of these processes is bounded from below. Consequently, any perturbation of an arbitrary wavelength cannot relax faster than an exponential function $e^{-{\mu }_{l}t}$, where ${\mu }_l$ is the lower bound.

This article is organized as follows. In Section 2 and Section 3, we discuss one-velocity models. These are processes characterized by a single velocity in norm, where only the direction of motion may change. These processes may seem, at a first glance, to be a simplistic physical model. However, in reality, they represent an important physical situation related to photon propagation in material media. In fact in this scenario, the transition in the velocity direction of the photon is solely dictated by scattering events [35].

In Section 4, we extend the analysis to processes characterized by a continuous structure of stochastic states. While the main qualitative result obtained in Section 2 and Section 3 is also confirmed in this case, novel and qualitatively interesting features arise, such as the incompleteness of the spectrum, and the fact that the eigenvalue spectrum can be defined solely for particular wave vectors (in the case discussed, these belong to an interval). The consequence of the latter property is that the relaxation decay in density dynamics is controlled not only by the point spectrum but also by the essential component of the spectrum. For a mathematical definition of the point and the essential spectra, the reader is referred to the classical monograph by Kato [36]. Throughout this article, we adopt an “operational” interpretation of the essential spectrum as the spectral complement to the point spectrum: an eigenvalue of the linear operator A defined in a functional space F (say the space of square summable functions) belongs to the essential spectrum if the corresponding eigenfunction does not belong to the functional space F.

Given the strong analogy between Poisson–Kac processes and quantum mechanics as provided by the seminal paper by Gaveau et al. [37], we explore an extension of the Gaveau model based on the functional structure of the model developed in Section 4, which possesses a continuum of velocity states, focusing on its spectral properties. This provides an interesting and exactly solvable example of a prototypical sub-quantum theory, in line with David Bohm [38,39], whose theory’s emergent properties coincide with the classical Schrödinger equation.

2. One-Velocity GPK Processes

Consider the simplest GPK process on the real line possessing two states $s=\{1,2\}$, a uniform transition rate ${\lambda }_0$ and a transition probability matrix $\mathbf{A}$ given by the following equation:

$$\mathbf{A}=\frac{1}{2}\left(\begin{array}{ccc}1+r& & 1+r\\ 1-r& & 1-r\end{array}\right),$$

(1)

where $r\in [0,1)$. The velocities in the two states are equal in absolute value and opposite in sign. We denote them as $b\left(1\right)=b(s=1)=b_0$, $b\left(2\right)=b(s=2)=-b_0$. This is therefore a one-velocity model. The stochastic model can be expressed as follows:

$$dx\left(t\right)=b\left({\chi }_2(t;{\lambda }_0,\mathbf{A})\right)dt,$$

(2)

where ${\chi }_2(t;{\lambda }_0,\mathbf{A})$ is a two-state finite Poisson process attaining the values $\{1,2\}$, i.e., a Markov process characterized by the transition rate ${\lambda }_0$ and by the transition probability matrix $\mathbf{A}$. The parameter r in Equation (1) determines a bias in the stochastic motion. The occurrence of two stochastic states implies that the statistical description of the process involves the partial probability densities $\mathbf{p}(x,t)=(p_1(x,t),p_2(x,t))$, and the SEO is given by $\partial \mathbf{p}(x,t)/\partial t=\mathcal{L}\left[\mathbf{p}\right(x,t\left)\right]$, where the infinitesimal generator $\mathcal{L}$ is [18]:

$$\mathcal{L}=\left(\begin{array}{ccc}-b_0\frac{\partial }{\partial x}-\tilde{\lambda }(1-r)& & \tilde{\lambda }(1+r)\\ & & \\ \tilde{\lambda }(1-r)& & b_0\frac{\partial }{\partial x}-\tilde{\lambda }(1+r)\end{array}\right),$$

(3)

while $\tilde{\lambda }={\lambda }_0/2$. Consider the associated eigenvalue problem $\mathcal{L}\left[\mathsf{\psi }\right(x\left)\right]=\mu \mathsf{\psi }\left(x\right)$. From the structure of $\mathcal{L}$ defined by Equation (3), the eigenfunctions are the imaginary exponentials $\mathsf{\psi }\left(x\right)=e^{ikx}({\psi }_1^0,{\psi }_2^0)$. Introducing the dimensionless quantities ${\mu }^{*}=\mu /\tilde{\lambda }$, $k^{*}=k{b}_0/\tilde{\lambda }$, the eigenvalues of $\mathcal{L}$ can be expressed by the following relation:

$${\mu }^{*}=-1+{\left[1-{\left(k^{*}\right)}^2-i2k^{*}r\right]}^{1/2},$$

(4)

where the square root should be interpreted in its multi-valued meaning, and $i=\sqrt{-1}$. The eigenvalue structure described by Equation (4) is depicted in Figure 1. More specifically: (i) the real part of the eigenvalues is lower bounded, i.e., $\mathrm{Re}\left[{\mu }^{*}\right]\ge -2$ (and of course, $\mathrm{Re}\left[{\mu }^{*}\right]\le 0$). For finite ${\lambda }_0$, this implies that the dynamics of any perturbation of the probability density $\delta \mathbf{p}(x,t)$, is lower bounded as $\left|\right|\delta \mathbf{p}(x,t)\left|\right|\ge C{e}^{-2t^{*}}$, where $\left|\right|·\left|\right|$ is any suitable norm (for instance, the $L^2$-norm) and $t^{*}=t\tilde{\lambda }$. For this process, therefore, ${\mu }_l={\lambda }_0$. (ii) For any wave vector k, two spectral branches exist. While in the symmetric case, $r=0$, the real parts of the spectral branches collapse into a single one for $k^{*}>1$, the effect of a bias, $r>0$, is to keep the two branches separate for any $k^{*}$. This phenomenon can be understood easily in the limit $k^{*}\to \infty$. Assuming $k^{*}\gg 1$ in Equation (4), we can write ${\mu }^{*}\simeq -1+{[-{\left(k^{*}\right)}^2-i2k^{*}r]}^{1/2}$. However, ${\left(k^{*}\right)}^2+2i{k}^{*}r={(k^{*}+ir)}^2+r^2\simeq {(k^{*}+ir)}^2$, and consequently, ${\mu }^{*}=-1\pm i(k^{*}+ir)$. This implies that the high wave vector region of the spectrum is well-approximated by the two eigenvalue branches ${\mu }^{*}=-(1+r)+i{k}^{*}$ and ${\mu }^{*}=-(1-r)-i{k}^{*}$, thus explaining the qualitative features depicted in Figure 1.

The latter result finds a direct interpretation in terms of the salient properties of the Green functions for this class of random motion. Consider the dynamics of an ensemble of particles moving according to Equation (2) and starting from the origin. Assume symmetric initial conditions, i.e., $p_1(x,0)=p_2(x,0)=\delta \left(x\right)/2$. Figure 2a depicts the evolution of the particle ensemble, expressed in terms of the overall probability density function, $p(x,t)=p_1(x,t)+p_2(x,t)$, obtained from stochastic simulations of $N_p=5\times {10}^7$ particles with $b_0=1$, ${\lambda }_0=1/2$, so that $D=b_0^2/2{\lambda }_0=1$. A characteristic feature of this class of finite-velocity stochastic motions is that the diagonal entries of the matrix-valued Green function are characterized by the superposition of a continuous, compactly supported component and of an impulsive contribution. Expressed in terms of the overall probability density $p(x,t)$ starting from an impulsive initial condition, this implies that $p(x,t)=p^c(x,t)+p^{\delta }(x,t)$, where $p^c(x,t)$ is the continuous term and $p^{\delta }(x,t)$ is the impulsive one, namely $p^{\delta }(x,t)=I_{+}\left(t\right)\delta (x-b_{0}t)+I_{-}\left(t\right)\delta (x+b_{0}t)$. The functions $I_{\pm }\left(t\right)$ progressively fade away as t increases [40,41,42]. The previous discussion suggests that the decay of the impulsive branches is controlled by the real part of the two asymptotic eigenvalues $\mathrm{Re}\left[{\mu }_1^{\infty }\right]=-\tilde{\lambda }(1+r)$ for $I_{-}\left(t\right)$ and $\mathrm{Re}\left[{\mu }_2^{\infty }\right]=-\tilde{\lambda }(1-r)$ for $I_{+}\left(t\right)$. This phenomenon is depicted in Figure 2b.

Next, we consider a two-dimensional GPK model possessing N stochastic states. This process can be characterized in terms of N partial probability densities $\mathbf{p}(\mathbf{x},t)={\left(p_{\alpha }(\mathbf{x},t)\right)}_{\alpha =1}^N$, satisfying the following general evolution equation [18]:

$$\frac{\partial p_{\alpha }}{\partial t}=-{\mathbf{b}}_{\alpha }·\nabla p_{\alpha }-{\lambda }_{\alpha }{p}_{\alpha }+\sum _{\beta =1}^N{A}_{\alpha ,\beta }{\lambda }_{\beta },p_{\beta }$$

(5)

with transition rates among states ${\lambda }_{\alpha }>0$, a probability transition matrix $A_{\alpha ,\beta }\ge 0$, where ${\sum }_{\alpha =1}^N{A}_{\alpha ,\beta }=1$ for $\alpha ,\beta =1,\cdots ,N$, and velocity vectors ${\mathbf{b}}_{\alpha }$. Moreover, in this case, we are interested in the spectral properties of the associated SEO (now vector-valued) acting on the partial densities defined by Equation (5). The eigenfunctions are planar waves, ${\psi }_{\alpha }\left(\mathbf{x}\right)=e^{i\mathbf{k}·\mathbf{x}}{\psi }_{\alpha }^0$, for $\alpha =1,\cdots ,N$. As regards the internal structure of the GPK model, we assume that all the transition rates are equal ${\lambda }_{\alpha }={\lambda }_0$, the transition probability matrix is uniform $A_{\alpha ,\beta }=1/N$, and the characteristic velocity vectors are ${\mathbf{b}}_{\alpha }=b_0{\mathsf{\beta }}_{\alpha }$, with ${\mathsf{\beta }}_{\alpha }=(cos{\phi }_{\alpha },sin{\phi }_{\alpha })$ and ${\phi }_{\alpha }=2\pi (\alpha -1)/N$ for $\alpha =1,\cdots ,N$. Introducing the non-dimensional quantities, ${\mathbf{k}}^{*}=\mathbf{k}{b}_0/{\lambda }_0$, ${\mu }^{*}=\mu /{\lambda }_0$, we can write the eigenvalue equation as follows:

$$\sum _{\beta =1}^N\left[A_{\alpha ,\beta }-(i{\mathbf{k}}^{*}·{\mathsf{\beta }}_{\alpha }+1){\delta }_{\alpha ,\beta }\right]{\psi }_{\beta }^0={\mu }^{*}{\psi }_{\alpha }^0,$$

(6)

where ${\delta }_{\alpha ,\beta }$ are the Kronecker symbols.

Figure 3 shows the real part of the eigenvalue spectrum of these models for two different values of N as a function of the norm $k^{*}=\left|{\mathbf{k}}^{*}\right|$ of the wave vector. The real part of the eigenvalue spectrum is lower bounded, specifically $\mathrm{Re}\left[{\mu }^{*}\right]\ge -1$. For this 2D GPK process, therefore, we obtain the same lower bound of the 1D biased model, ${\mu }_l={\lambda }_0$. The comparison of the data shown in the two panels suggests the existence of a well-defined spectral limit for $N\to \infty$, i.e., in the case where the velocities are uniformly distributed on the unit circumference. This limiting case, which corresponds to the Markov motions analyzed by Kolesnik [17], suggests that for a continuum of stochastic states, the spectral properties of the resulting SEO give rise to a much simpler structure for the spectral branches.

3. Lévy Walks

Owing to the analysis developed by Fedotov [29], subsequently elaborated in [30,31,32,33] and extended in a unitary theory of stochastic processes possessing bounded velocity in [34], a complete statistical description of an LW involves a system of partial probability densities, as for GPK processes. The only difference is that in the case of LW, we also need to parametrize the partial densities with respect to an internal parameter, the transition age. This is due to the more complex nature of the density function for the transition times that is no longer exponential (as for Poisson–Kac processes). We consider a one-dimensional LW switching between two states, corresponding to the directions of motion, keeping constant the absolute value of the velocity $b_0$. The statistical properties of such an LW are fully described by the partial probability densities $p_{\pm }(x,t;\tau )$. Setting the following equation:

$$\mathbf{p}(x,t;\tau )=\left(\begin{array}{c}{p}_{+}(x,t;\tau )\\ p_{-}(x,t;\tau )\end{array}\right),$$

(7)

the vector-valued density $\mathbf{p}(x,t;\tau )$ satisfies the evolution equation [29]

$$\frac{\partial \mathbf{p}}{\partial t}={\mathcal{M}}_{\tau }\left[\mathbf{p}\right]+{\mathcal{A}}_x\left[\mathbf{p}\right]={\mathcal{L}}_{\tau ,x}\left[\mathbf{p}\right],$$

(8)

where ${\mathcal{M}}_{\tau }$ represents the evolution operator for the transition age,

$${\mathcal{M}}_{\tau }=\left(\begin{array}{cc}-\frac{\partial }{\partial \tau }-\lambda \left(\tau \right)& 0\\ 0& -\frac{\partial }{\partial \tau }-\lambda \left(\tau \right)\end{array}\right)=-{\mathbf{I}}_2\left(\frac{\partial }{\partial \tau }+\lambda \left(\tau \right)\right),$$

(9)

${\mathbf{I}}_2$ is the $2\times 2$ identity matrix, and ${\mathcal{A}}_x$ represents the advection operator:

$${\mathcal{A}}_x=\left(\begin{array}{cc}-b_0\frac{\partial }{\partial x}& 0\\ 0& b_0\frac{\partial }{\partial x}\end{array}\right)=-{\mathsf{\sigma }}_z{b}_0\frac{\partial }{\partial x},$$

(10)

with ${\mathsf{\sigma }}_z=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)$, a Pauli matrix. Both operators act on the two partial probability density waves $p_{\pm }$. Assuming that particles invert their direction of motion at each transition, Equation (8) is equipped with the following boundary condition:

$$\mathbf{p}(x,t,0)={\sigma }_x{\int }_0^{\infty }\lambda \left(\tau \right)\mathbf{p}(x,t,\tau )d\tau ,$$

(11)

where ${\mathsf{\sigma }}_x$ is the Pauli matrix ${\mathsf{\sigma }}_x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$. We now consider the spectral structure underlying the SEO, that is, we solve the eigenvalue problem:

$${\mathcal{L}}_{\tau ,x}\left[\mathsf{\psi }(x,\tau )\right]=\mu \mathsf{\psi }(x,\tau ),\mathsf{\psi }(x,\tau )=\left(\begin{array}{c}{\psi }_{+}(x,\tau )\\ {\psi }_{-}(x,\tau )\end{array}\right).$$

(12)

The eigenfunction $\mathsf{\psi }(x,\tau )$ can be expressed as planar waves with respect to the spatial coordinate, i.e., $\mathsf{\psi }(x,\tau )=e^{ikx}\mathsf{\varphi }\left(\tau \right)$, while the dependence on the transition time $\tau$ is accounted for by the vector-valued function $\mathsf{\varphi }\left(\tau \right)=({\varphi }_{+}\left(\tau \right),{\varphi }_{-}\left(\tau \right))$. Its two components satisfy the following equations:

$$\frac{d{\varphi }_{\pm }\left(\tau \right)}{d\tau }=-\left(\mu \pm ik{b}_0+\lambda \left(\tau \right)\right){\varphi }_{\pm }\left(\tau \right),$$

(13)

and thus,

$${\varphi }_{\pm }\left(\tau \right)=A_{\pm }exp\left[-\mu \tau \mp ik{b}_0\tau -{\int }_0^{\tau }\lambda \left(\theta \right)d\theta \right],$$

(14)

where $A_{\pm }$ are two complex-valued constants.

Imposing the boundary conditions (11) on the function $\mathsf{\varphi }$, we find that the two constants $A_{\pm }$ satisfy the following relation:

$$A_{\pm }=A_{\mp }{\int }_0^{\infty }T\left(\tau \right)e^{-(\mu \pm ik{b}_0)\tau }d\tau ,$$

(15)

where $T\left(\tau \right)=\lambda \left(\tau \right)exp\left[-{\int }_0^{\tau }\lambda \left(\theta \right)d\theta \right]$ is the probability density function for the transition times. Eliminating $A_{\pm }$ from the two in Equation (15), we obtain the characteristic equation for the eigenvalues of ${\mathcal{L}}_{\tau ,x}$, given as follows:

$$\left[{\int }_0^{\infty }T\left(\tau \right)e^{-(\mu +ik{b}_0)\tau }d\tau \right]\left[{\int }_0^{\infty }T\left({\tau }^{\prime }\right)e^{-(\mu -ik{b}_0){\tau }^{\prime }}d{\tau }^{\prime }\right]=1.$$

(16)

For $k=0$, one recovers the spectral properties of the “renewal mechanism” of the LW. In this case, Equation (16) reduces to the following:

$${\left[{\int }_0^{\infty }T\left(\tau \right)e^{-\mu \tau }d\tau \right]}^2=1.$$

(17)

Equation (17) can be interpreted as follows: the spectral structure of the renewal equation of an LW corresponds to the zeroes of the equations $\widehat{T}\left(\mu \right)=\pm 1$, where $\widehat{T}\left(\mu \right)$ represents the Laplace transform (of argument $\mu$) of the transition-time probability density $T\left(\tau \right)$. For real $\mu$, only the positive determination of the latter equation makes sense, and it can be clearly noted that there exists a unique real solution for $\widehat{T}\left(\mu \right)=1$, namely $\mu =0$, that corresponds to the Frobenius (conservation) eigenvalue of the process.

In the rest of this section, we consider three typical LW models and solve for their spectral properties.

3.1. Gamma-Distributed Transition Times

For this model, the transition-time density is expressed by the following equation:

$$T\left(\tau \right)=\frac{{\beta }^{\alpha }}{\mathsf{\Gamma }\left(\alpha \right)}{\tau }^{\alpha -1}{e}^{-\beta \tau },$$

(18)

with $\alpha \in (0,\infty )$, $\beta >0$. Since

$$\widehat{T}\left(s\right)=\frac{{\beta }^{\alpha }}{{(s+\beta )}^{\alpha }},$$

(19)

it follows that Equation (16) reduces to

$${\left[{(\mu +\beta )}^2+k^2{b}_0^2\right]}^{\alpha }={\beta }^{2\alpha },$$

(20)

or equivalently,

$${(\mu +\beta )}^2+k^2{b}_0^2={\beta }^2{e}^{i{\phi }_{\alpha ,h}},$$

(21)

where ${\phi }_{\alpha ,h}=2\pi h/\alpha$, and h is an integer, $h=0,1,\cdots$. Two cases may occur: (i) if $\alpha =P/Q$ is rational, with $P,Q$ both being integers, there are Q distinct values for $e^{i{\phi }_{\alpha ,h}}$ that correspond to $h=0,1,\cdots ,Q$; (ii) if $\alpha$ is irrational, there is a countable number of distinct values for $e^{i{\phi }_{\alpha ,h}}$, each corresponding to a different spectral branch.

The spectral properties of the renewal mechanism (corresponding to $k=0$) follow from Equation (21):

$${\mu }_h=-\beta \pm \beta e^{i{\phi }_{\alpha ,h}/2},$$

(22)

while the eigenvalues are generally complex-valued. In particular, if $\alpha$ is irrational, $({\mu }_h+\beta )/\beta$, $h=0,1,\cdots$, densely fill the unit circumference.

Next, we consider the case $k\ne 0$. Solving the quadratic Equation (21), one obtains the following:

$${\mu }_h=-\beta \pm {\left[{\beta }^2{e}^{i{\phi }_{\alpha ,h}}-k^2{b}_0^2\right]}^{1/2}=-\beta \pm \sqrt{z_h},$$

(23)

where $z_h=\left[{\beta }^{2}cos\left({\phi }_{\alpha ,h}\right)-k^2{b}_0^2\right]+i{\beta }^{2}sin\left({\phi }_{\alpha ,h}\right)={\rho }_h{e}^{i{\omega }_h}$. Since the two determinations of $\sqrt{z_h}$ are $\pm {\rho }_h^{1/2}{e}^{i{\omega }_h/2}$, it follows that:

$${\mu }_h=-\beta \pm {\rho }_h^{1/2}\left[cos\left(\frac{{\omega }_h}{2}\right)+isin\left(\frac{{\omega }_h}{2}\right)\right].$$

(24)

It can be clearly shown that the real part of the spectrum is lower bounded. Moreover, the explicit calculation of the eigenvalues indicates that:

$$-2\beta \le \mathrm{Re}\left[{\mu }_h\right]\le 0.$$

(25)

In order to give some examples, Figure 4 and Figure 5 depict the real and imaginary part of the spectrum in two typical cases: $\alpha =3$ (integer) and $\alpha =\pi$ (irrational). The other parameters are set to $\beta =b_0=1$. In the irrational case, solely the the spectral branches corresponding to $h=0,\cdots ,19$, are depicted.

3.2. Superdiffusive Lévy Walk

For LWs defined by a transition-time probability density possessing a long-term power-law scaling, such as those associated with the distribution function $F_{\tau }\left(\tau \right)={\int }_0^{\tau }T\left({\tau }^{\prime }\right)d{\tau }^{\prime }=1-E_{\nu }(-{(\tau /{\tau }_0)}^{\nu })$, where $E_{\nu }\left(z\right)={\sum }_{k=0}^{\infty }{z}^k/\mathsf{\Gamma }(\nu k+1)$ is the Mittag–Leffler function, $0<\nu <1$, ${\tau }_0>0$, and $\mathsf{\Gamma }\left(z\right)$ the Gamma function [43,44], the Laplace transform of $T\left(\tau \right)$ is thus given by the following equation:

$$\widehat{T}\left(s\right)=\frac{1}{1+{\left({\tau }_{0}s\right)}^{\nu }},$$

(26)

and Equation (16) provides:

$$\left[1+{\tau }_0^{\nu }{(\mu +ik{b}_0)}^{\nu }\right]\left[1+{\tau }_0^{\nu }{(\mu -ik{b}_0)}^{\nu }\right]=1.$$

(27)

To solve this equation, we set the following:

$$\mu +ik{b}_0=\mathrm{Re}\left[\mu \right]+i\left(\mathrm{Im}\left[\mu \right]\pm k{b}_0\right)=\rho e^{i{\varphi }_{\pm }},$$

(28)

where $\rho \ge \mathrm{Re}\left[\mu \right]$. Equation (27) can be thus rewritten as follows:

$$1+{\tau }_0^{\nu }{\rho }^{\nu }\left(e^{i\nu {\varphi }_{+}}+e^{i\nu {\varphi }_{-}}\right)+{\tau }_0^{2\nu }{\rho }^{2\nu }{e}^{i\nu ({\varphi }_{+}+{\varphi }_{-})}=1,$$

(29)

which can be further reduced to the equation below:

$${\tau }_0^{\nu }{\rho }^{\nu }+(e^{-i\nu {\varphi }_{+}}+e^{-i\nu {\varphi }_{-}})=0.$$

(30)

This implies the following relation:

$$\left|\mathrm{Re}\left[\mu \right]\right|\le \rho \le \frac{2^{1/\nu }}{{\tau }_0},$$

(31)

which proves the boundedness of the real part of the spectrum.

3.3. Power-Law Distributed Transition Times

We consider as a third example the classical case of an LW (8)–(11) defined by the the transition-time probability density $T\left(\tau \right)=\xi /{(1+\tau )}^{\xi +1}$, where $\xi >0$. Its Laplace transform is given by the following equation:

$$\widehat{T}\left(\mu \right)=\xi {\int }_0^{\infty }\frac{e^{-\mu \tau }}{{(1+\tau )}^{\xi +1}}d\tau =\xi e^{\mu }{Ei}_{\xi +1}\left(\mu \right),$$

(32)

where ${\mathrm{Ei}}_{\xi +1}\left(z\right)$ is the Exponential Integral function of order $\xi +1$,

$${\mathrm{Ei}}_{\nu }\left(z\right)={\int }_1^{\infty }\frac{e^{-zy}}{y^{\nu }}dy,$$

(33)

In this case, it is not possible to derive the structure of the eigenvalue branches in closed form. Nevertheless, we can still obtain spectral bounds using asymptotic analysis.

Let us prove ad absurdum that there exists a lower bound for the real part of the eigenvalue spectrum. Let us suppose the opposite, namely that there are eigenvalues with arbitrarily large real parts, which implies arbitrarily large values for $\left|\mu \right|$. In this case, we can use Equation (32), and the asymptotic expansion for the Exponential Integral function ${\mathrm{Ei}}_{\nu }\left(z\right)$ for complex z [45]:

$${\mathrm{Ei}}_{\nu }\left(z\right)=\frac{e^{-z}}{z}\left[1-\frac{\nu }{z}+O\left(\frac{1}{z^2}\right)\right],$$

(34)

which implies the following:

$$\widehat{T}(\mu \pm ik{b}_0)=\frac{\xi }{\mu \pm ik{b}_0}+{\left.O\left(\frac{1}{s^3}\right)\right|}_{s=\mu \pm ik{b}_0}.$$

(35)

For large $\mathrm{Re}\left[\mu \right]$, one obtains the following from Equation (16):

$$\frac{{\xi }^2}{(\mu +ik{b}_0)(\mu -ik{b}_0)}=1+O\left(\frac{1}{\mathrm{Re}{\left[\mu \right]}^4}\right).$$

(36)

By hypothesis, the real part of the eigenvalues can attain arbitrarily large values so that the higher-order terms in Equation (36) can be neglected, reducing it to the equation below:

$${\mu }^2+k^2{b}_0^2-{\xi }^2=0.$$

(37)

The solutions of this equation are as follows:

$$\mu =\left\{\begin{array}{ccc}\pm \sqrt{{\xi }^2-k^2{b}_0^2}\hfill & \hfill & \left|k\right|<\xi /b_0\hfill \\ \pm i\sqrt{k^2{b}_0^2-{\xi }^2}\hfill & \hfill & \left|k\right|>\xi /b_0.\hfill \end{array}\right.$$

(38)

It thus follows from Equation (38) that the real part of the spectrum is lower-bounded by $-\xi /b_0$, contradicting the original hypothesis. The proof is completed.

From the properties of the Laplace transforms, namely ${lim}_{\mathrm{Re}\left[\mu \right]\to \infty }\widehat{T}\left(\mu \right)=0$, and applying the same analysis developed above, we can show that all the one-dimensional LWs, whose spectrum fulfills Equation (16), are characterized by a lower bound $-{\mu }_l$, with ${\mu }_l>0$ for $\mathrm{Re}\left[\mu \right]$.

4. GPK Process with a Continuum of States

A common feature of one-velocity models is the decomposition of the Green function into a continuous part and an impulsive part, as discussed in [40,42] and illustrated in Figure 2.

A question naturally arises, as to whether it would be possible to define a stochastic process with a finite propagation velocity, for which the corresponding matrix-valued Green function solely admits the continuous part but not an impulsive contribution, similarly to what happens to the solution of the parabolic diffusion equation. While the answer to this question is negative (for reasons that are addressed below), it is indeed possible to construct processes such that if the overall initial density is impulsive, say $p(x,0)=\delta \left(x\right)$, then for any $t>0$, the overall density function $p(x,t)$ is smooth and compactly supported, i.e., it does not not possess any impulsive term. This construction is interesting when further addressing the spectral properties of these models, as the sudden annihilation of the impulsive initial conditions necessarily implies a spectral counterpart in the high-wave vector limit. Examples of these models have been discussed in [46].

Processes that possess this property implies a continuous set of stochastic states possessing velocities defined in some bounded and continuous set (in the present case, an interval) and smooth transitions among them. In this case, the velocity itself can be used to parametrize the stochastic states. A typical example is given by the stochastic dynamics:

$$dx\left(t\right)=b_0\Xi (t;\lambda ;K)dt,$$

(39)

where $\Xi (t;\lambda ,K)$ is a Poisson field defined in [47], and attaining values in $[-1,1]$, i.e., a stochastic process such that the densities $\widehat{P}(\beta ,t)$ of the probabilities, where $\widehat{P}(\beta ,t)d\beta =\mathrm{Prob}[\Xi \left(t\right)\in (\beta ,\beta +d\beta )]$, satisfy a Markov dynamics:

$$\frac{\partial \widehat{P}(\beta ,t)}{\partial t}=-\lambda \left(\beta \right)\widehat{P}(\beta ,t)+{\int }_{-1}^{1}K(\beta ,{\beta }^{\prime })\lambda \left({\beta }^{\prime }\right)\widehat{P}({\beta }^{\prime },t)d{\beta }^{\prime },$$

(40)

where $\lambda \left(\beta \right)$ is a smooth positive function expressing the transition rate from the stochastic state $\beta$, and $K(\beta ,{\beta }^{\prime })$ a smooth stochastic kernel, with $K(\beta ,{\beta }^{\prime })\ge 0$, and

$${\int }_{-1}^{1}K(\beta ,{\beta }^{\prime })d\beta =1,$$

(41)

for all ${\beta }^{\prime }\in [-1,1]$. The latter accounts for the transition probabilities from state ${\beta }^{\prime }$ to all the other states. Henceforth, we consider the simplest case of uniform $\lambda \left(\beta \right)$ and $K(\beta ,{\beta }^{\prime })$, i.e.,:

$$\lambda \left(\beta \right)={\lambda }_0,K(\beta ,{\beta }^{\prime })=\frac{1}{2}.$$

(42)

The statistical description of the stochastic motion defined by Equation (39) involves the family of partial probability densities $p(x,t;\beta )$, continuously parametrized with respect to $\beta \in [-1,1]$ and satisfying the following equation:

$$\frac{\partial p(x,t;\beta )}{\partial t}=-b_0\beta \frac{\partial p(x,t;\beta )}{\partial x}-{\lambda }_{0}p(x,t;\beta )+\frac{{\lambda }_0}{2}{\int }_{-1}^{1}p(x,t;{\beta }^{\prime })d{\beta }^{\prime }.$$

(43)

In this model, the overall probability density $P(x,t)$ for the particle position x at time t and the associated flux $J(x,t)$ are given by the equation below:

$$P(x,t)={\int }_{-1}^{1}p(x,t;\beta )d\beta ,J(x,t)={\int }_{-1}^1\beta p(x,t;\beta )d\beta .$$

(44)

Using the method outlined in [18,47,48], it follows that the Kac limit of Equation (43), i.e., the limit for the unbounded propagation velocity and the transition rate, $b_0,{\lambda }_0\to \infty$, keeping fixed the nominal diffusivity $D_{\mathrm{nom}}=b_0^2/2{\lambda }_0$ is given by the diffusion equation for $P(x,t)$:

$$\frac{\partial P(x,t)}{\partial t}=D\frac{{\partial }^{2}P(x,t)}{\partial x^2},$$

(45)

with a value of the effective diffusivity D equal to $D=2D_{\mathrm{nom}}/3$. Moreover, for any finite value of $b_0$ and ${\lambda }_0$, the long-term solutions of Equation (43) approach those of the parabolic Equation (45).

Consider the solutions of Equation (43), starting from a spatially impulsive initial distribution $P(x,0)=\delta \left(x\right)$. This condition does not completely specify the initial state of the system, as the initial preparation with respect to the internal parametrization $\beta$ should also be defined, i.e., the whole structure of the partial densities $p(x,0;\beta )={\tilde{p}}_0\left(\beta \right)\delta \left(x\right)$, $\beta \in [-1:1]$ at $t=0$ should be specified. Two situations can occur: (i) if ${\tilde{p}}_0\left(\beta \right)=\delta (\beta -{\beta }^{*})$ admits an impulsive component at any $\beta ={\beta }^{*}$, then $P(x,t)$ for $t>0$ is characterized by an impulsive contribution centered at $x={\beta }^{*}t$, propagating with velocity ${\beta }^{*}$ superimposed onto a continuous distribution that derives from the recombination mechanism among the partial waves. Conversely, (ii) if ${\tilde{p}}_0\left(\beta \right)$ is smooth, i.e., no impulsive initial term is present, then $P(x,t)$ for $t>0$ is also a smooth function of x for any $t>0$. The phenomena outlined above are depicted in Figure 6a,b. Data have been obtained from stochastic simulations using an ensemble of $N_p=5\times {10}^7$ walkers initially placed at $x=0$. In the case of the data depicted in panel (a), the initial distribution with respect to the internal state variable is impulsive and centered at ${\beta }^{*}=1/2$. Conversely, the initial condition for the data depicted in panel (b) is uniform over $\beta$, i.e., ${\tilde{p}}_0\left(\beta \right)=1/2$, $\beta \in [-1,1]$.

This example clearly indicates that even in the presence of a continuum of internal states, any initial preparation that is impulsive with respect to the internal-state parametrization $\beta$ determines an overall density $P(x,t)$ that contains an impulsive contribution with respect to x, propagating at constant speed. It is interesting to analyze how this phenomenology can be interpreted in the light of the spectral properties of the associated evolution operator.

4.1. Spectral Properties

Consider the spectral properties of the infinitesimal generator of the density dynamics (43):

$${\mathcal{L}}_{\beta }\left[\varphi (x,\beta )\right]=-b_0\beta \frac{\partial \varphi (x,\beta )}{\partial x}-{\lambda }_0\varphi (x,\beta )+\frac{{\lambda }_0}{2}{\int }_{-1}^1\varphi (x,{\beta }^{\prime })d{\beta }^{\prime }.$$

(46)

The eigenvalue equation for ${\mathcal{L}}_{\beta }$ is as follows:

$${\mathcal{L}}_{\beta }\left[\psi (x,\beta )\right]=\mu \psi (x,\beta ),$$

(47)

which admits the eigenfunctions of the following form:

$$\psi (x,\beta )=e^{ikx}{\psi }_k\left(\beta \right).$$

(48)

Introducing the dimensionless quantities ${\mu }^{*}=\mu /{\lambda }_0$, $k^{*}=k{b}_0/{\lambda }_0$, the substitution of Equation (48) into Equations (46) and (47) provides the following equation:

$${\psi }_k\left(\beta \right)=\frac{1}{2(1+{\mu }^{*}+i{k}^{*}\beta )}{\int }_{-1}^1{\psi }_k\left({\beta }^{\prime }\right)d{\beta }^{\prime }.$$

(49)

Consequently, integrating over $\beta$ and assuming ${\int }_{-1}^1{\psi }_k\left(\beta \right)d\beta \ne 0$, the characteristic equation for the eigenvalues is as follows:

$${\int }_{-1}^1\frac{d\beta }{1+{\mu }^{*}+i{k}^{*}\beta }=2.$$

(50)

Setting ${\mu }^{*}=r+i\omega$ and expliciting the integral in Equation (50), one obtains the following equation:

$$\arctan \left(\frac{\omega +k^{*}}{1+r}\right)-\arctan \left(\frac{\omega -k^{*}}{1+r}\right)-\frac{i}{2}\log \left[\frac{{(1+r)}^2+{(\omega +k^{*})}^2}{{(1+r)}^2+{(\omega -k^{*})}^2}\right]=2k^{*}.$$

(51)

Imposing the imaginary part to be null, one obtains $\omega =0$, i.e., the eigenvalues are purely real, ${\mu }^{*}=r$. Enforcing this property back in Equation (51), a simpler equation for the real part r follows:

$$k^{*}=\arctan \left(\frac{k^{*}}{1+r}\right).$$

(52)

The solution of Equation (52) can be easily obtained by introducing the auxiliary variable $z=k^{*}/(1+r)$, so that Equation (52) becomes $(1+r)z=\arctan \left(z\right)$. Varying $z\in [0,\infty )$, one gets the continuous branch of eigenvalues:

$$r=-1+\frac{\arctan \left(z\right)}{z}.$$

(53)

Given z, the corresponding wavenumber $k^{*}$ follows from the definition of z, namely $k^{*}=(1+r)z$. Before addressing the properties of this spectral branch, let us show that there is no other eigenvalue. To show this, we have to consider the only remaining condition:

$${\int }_{-1}^1{\psi }_k\left(\beta \right)d\beta =0.$$

(54)

In this case, the eigenvalue problem reduces to $({\mu }^{*}+1+i{k}^{*}\beta ){\psi }_k\left(\beta \right)=0$, which admits no solution for a ${\mu }^{*}$ independent of $\beta$. Consequently, the spectrum reduces to the spectral branch defined by Equations (52) and (53). This branch is defined for $|k^{*}|<k_c=\pi /2$, as depicted in Figure 7. In this interval of wave vectors, the associated eigenfunctions defined by Equation (49) are complex-valued and can be expressed as follows:

$${\psi }_k\left(\beta \right)=A\left[{\varphi }_k\left(\beta \right)+i{\varphi }_k(-\beta )\right],$$

(55)

where A, is a complex-valued constant and ${\varphi }_k\left(\beta \right)$ is a real-valued function defined as:

$${\varphi }_k\left(\beta \right)=C\frac{(1+{\mu }^{*})+k^{*}\beta }{{(1+{\mu }^{*})}^2+{\left(k^{*}\beta \right)}^2},$$

(56)

where the normalization constant C is such that $\left|\right|{\varphi }_k\left(\beta \right){\left|\right|}_{L^2}={\int }_{-1}^1{\varphi }_k^2\left(\beta \right)d\beta =1$. Expression (55) with ${\varphi }_k\left(\beta \right)$ given by Equation (56) follows from Equation (49), i.e., ${\psi }_k\left(\beta \right)=B/[(1+{\mu }^{*})+i{k}^{*}\beta ]$, setting $B=1+i$, modulo being an arbitrary multiplicative constant. Figure 8 depicts the shape of ${\varphi }_k\left(\beta \right)$ at $k^{*}=1.5$, close to the break-up point $k_c\simeq 1.57$.

From the functional form of the eigenfunctions (55) and (56), the dynamical mechanism associated with the spectral break-up becomes clear. As the eigenvalue ${\mu }^{*}$ approaches $-1$ from below, the associated eigenfunction develops a singularity at $\beta =0$, ${\psi }_k\left(\beta \right)\to 1/\beta$, and consequently the critical point $k_c$ corresponds to the bifurcation point where ${\psi }_k\left(\beta \right)$ ceases to be summable.

The occurrence of a single eigenfunction, restricted solely to the interval $[-k_c,k_c]$, is manifestly incomplete in order to represent a generic (square-summable) function of the internal variable $\beta$ in the interval $[-1,1]$. It is therefore interesting to analyze the relaxation properties of the operator ${\mathcal{L}}_{\beta }$, restricted to a planar wave mode $p(x,t;\beta )=e^{ikx}\zeta (\beta ,t)$, i.e., to consider the following dynamics:

$$\frac{\partial \zeta (\beta ,t)}{\partial t}=-i{k}^{*}\beta \zeta (\beta ,t)-\zeta (\beta ,t)+\frac{1}{2}{\int }_{-1}^1\zeta ({\beta }^{\prime },t)d{\beta }^{\prime },$$

(57)

where $\zeta (\beta ,t)$ is complex-valued, and the time variable t has been made non-dimensional by normalizing it with ${\lambda }_0$. It thus follows from Equation (57) that:

$$\begin{array}{ccc}\hfill \frac{{\partial \left|\zeta (\beta ,t)\right|}^2}{\partial t}& =& \overline{\zeta }(\beta ,t)\frac{\partial \zeta (\beta ,t)}{\partial t}+\zeta (\beta ,t)\frac{\partial \overline{\zeta }(\beta ,t)}{\partial t}=-2{\left|\zeta (\beta ,t)\right|}^2\hfill \\ & +& \mathrm{Re}\left[\overline{\zeta }(\beta ,t){\int }_{-1}^1\zeta ({\beta }^{\prime },t)d{\beta }^{\prime }\right],\hfill \end{array}$$

(58)

where $\mathrm{Re}[·]$ is the real part of the argument, and $\overline{\zeta }(\beta ,t)$ the complex conjugate of $\zeta (\beta ,t)$. Integrating over $\beta$, the dynamics of the $L^2$-norm of $\zeta (\beta ,t)$, ${\left|\right|\zeta \left|\right|}^2\left(t\right)={\int }_{-1}^1{\left|\zeta (\beta ,t)\right|}^{2}d\beta$ is as follows:

$$\frac{{d\left|\right|\zeta \left|\right|}^2\left(t\right)}{dt}=-{2\left|\right|\zeta \left|\right|}^2\left(t\right)+{\left|Z\left(t\right)\right|}^2,$$

(59)

where $Z\left(t\right)={\int }_{-1}^1\zeta (\beta ,t)d\beta$. As ${\left|Z\left(t\right)\right|}^2\ge 0$, it follows from Equation (59) that:

$$\frac{{d\left|\right|\zeta \left|\right|}^2\left(t\right)}{dt}\ge -{2\left|\right|\zeta \left|\right|}^2\left(t\right)\Rightarrow \left|\right|\zeta \left|\right|\left(t\right)\ge \left|\right|\zeta \left|\right|\left(0\right)e^{-t}.$$

(60)

Equation (60) indicates that all the planar-wave modes cannot decay to zero faster than $e^{-{\gamma }^{*}t}$, with ${\gamma }^{*}=1$ representing the maximum relaxation rate. This result is valid for any $k^{*}$, below and above the critical value $k_c$.

Consequently, in the statistical evolution of the process, the point spectrum plays a marginal role, especially for high wave vectors, as the dynamics is controlled by the essential spectral component. From Equation (60), it also follows that the essential spectrum is lower bounded, and this fact is expressed by the property ${\gamma }^{*}=1$. The typical evolution of a high-wave vector excitation above the critical value $k_c$ is depicted in Figure 9, showing the progressive development of a complex structure in $\beta$, which does not converge to any proper eigenfunction.

By expanding $\zeta (\beta ,t)$ in the truncated Fourier series with respect to $\beta$, $\zeta (\beta ,t)=$${\sum }_{n=-N}^N{\zeta }_n\left(t\right)e^{in\pi \beta }$, keeping N sufficiently large, the eigenvalue spectrum can be calculated by solving the eigenvalue problem ${\sum }_{k=-N}^N{M}_{n,k}{\zeta }_k={\mu }^{*}{\zeta }_n$, where $M_{0,0}=0$, $M_{n,k}=-k^{*}d(k-n)/2-{\delta }_{n,k}$, for $\left|n\right|+\left|k\right|>0$, where $d\left(k\right)=2{(-1)}^k/\left(k\pi \right)$, $\left|k\right|>0$. The direct numerical calculation of the spectrum, obtained numerically by setting $N=5000$, provides for $k^{*}<k_c$, only a single eigenvalue in the point spectrum, as already found. Meanwhile, for all the values of $k^{*}$, the essential spectrum is formed by essential eigenvalues ${\mu }_{\mathrm{ess}}^{*}$ possessing a real part that is identical and equal to $-1$, and an imaginary part distributed in an uniform way within the interval $[-k^{*},k^{*}]$. This phenomenon is depicted in Figure 10.

The overall dynamics of the process involves both the spatial evolution (i.e., the dynamics with respect to the variable x), parametrized by the wave vector k (upon a Fourier transform), and the internal dynamics with respect to the variable $\beta$, describing the recombination process among the partial density waves. Below the critical threshold $k_c$, the unique point spectrum eigenvalue controls the long-term exponential relaxation of the solutions of Equation (57). However, the initial/intermediate decay is sensitive to the relaxation exponent ${\gamma }^{*}$ that pertains to the essential spectrum. This initial/intermediate behavior becomes more evident, starting from initial conditions $\zeta (\beta ,0)$ that possess high-frequency components with respect to the internal variable $\beta$. This phenomenon can be clearly observed from the data in Figure 11, which shows the decay of the $L^2$-norm $\left|\right|\zeta \left|\right|\left(t\right)$ of $\zeta (\beta ,t)$ at $k^{*}=1$, starting from differential initial sinusoidal profiles, $\zeta (\beta ,0)=sin\left(\nu \pi \beta \right)$, $\nu$ integer.

For $\nu \sim O\left(1\right)$, only the long-term exponential relaxation can be observed (line a), the exponent of which is the single point-spectrum exponent $r\simeq 0.356$. As $\nu$ increases, $\nu >{10}^1$, a crossover between the initial/intermediate exponential scaling that depends on the essential component of the spectrum $\left|\right|\zeta \left|\right|\left(t\right)\sim e^{-{\gamma }^{*}t}$ and on the long-term relaxation $\left|\right|\zeta \left|\right|\left(t\right)\sim e^{-rt}$ occurs.

As a final observation, let us consider the decay of the impulsive component of the overall density function for an initial preparation of the system, with particles possessing one and the same initial velocity, i.e., $p(x,t;\beta )=\delta \left(x\right)\delta (\beta -{\beta }^{*})$. This situation corresponds qualitatively to the evolution of the density profiles depicted in Figure 2a for the simpler one-velocity model. In this case, the overall density $P(x,t)$ at time t is the superposition of a continuous distribution $P_c(x,t)$ and of an impulsive component centered at $x=b_0{\beta }^{*}t$:

$$P(x,t;)=P_c(x,t)+I_{{\beta }^{*}}\left(t\right)\delta (x-b_0{\beta }^{*}t).$$

(61)

The impulsive component $I_{{\beta }^{*}}\left(t\right)$, translating in space at velocity ${\beta }^{*}$, corresponds to the fraction of particles that remain in the initial state ${\beta }^{*}$, i.e., not experiencing any transition up to time t. Because of its physical meaning, it is clear that the relaxation of $I_{{\beta }^{*}}\left(t\right)$ should not depend on the spatial component of the dynamics of $p(x,t;\beta )$, i.e., should be independent of the wave vector k. Furthermore, it follows from the above observation that this is exclusively a property of the Markov recombination mechanism of the partial density waves, as described by the internal state dynamics of Equation (40), and associated with the internal state operator $\widehat{\mathcal{L}}$:

$$\widehat{\mathcal{L}}\left[f\left(\beta \right)\right]=-{\lambda }_{0}f\left(\beta \right)+\frac{{\lambda }_0}{2}{\int }_{-1}^{1}f\left({\beta }^{\prime }\right)d{\beta }^{\prime }.$$

(62)

The spectrum of $\widehat{\mathcal{L}}$ is completely degenerate: it possesses only the single eigenvalue $\mu ={\lambda }_0$ with countable multiplicity, as any function $\pi \left(\beta \right)$ possessing a zero mean in $[-1,1]$ is an eigenfunction associated with it. More specifically, a basis for this eigenspace is given by the sinusoidal functions ${\pi }_{\nu }\left(\beta \right)=sin\left(\nu \pi \beta \right)$, where $\nu =1,2,\dots$ This eigenproperty of the recombination mechanism, when embedded into the spatio-temporal evolution of the partial density waves $p(x,t,\beta )$ defined by Equation (43), produces results incompatible with the translation operator $-b_0\beta \nabla$ in order to determine any proper eigenstructure, thus contributing to the essential part of the spectrum of the operator ${\mathcal{L}}_{\beta }$ associated with the spatial-temporal evolution of the partial density waves.

It follows from the above discussion that the impulsive part of the dynamics $I_{{\beta }^{*}}\left(t\right)$ entering Equation (62) should decay exponentially as follows:

$$I_{{\beta }^{*}}\left(t\right)\sim e^{-{\lambda }_{0}t}.$$

(63)

corresponding to the decay defined by the essential spectral component. This phenomenon is depicted in Figure 12. The numerical data have been obtained from the stochastic simulation of an ensemble of $N_p={10}^8$ particles evolving in space according to Equation (43), initially placed at $x=0$ and all possessing the same velocity $b_0{\beta }^{*}$. Therefore, this preparation of the system corresponds to the initial condition for the partial density waves expressed by $p(x,0;\beta )=\delta \left(x\right)\delta (\beta -{\beta }^{*})$. By tracking the evolution of the particle system, the intensity of the impulsive peak $I_{{\beta }^{*}}\left(t\right)$ can be determined.

Data in Figure 12 refer to different values of $b_0$, keeping them fixed at the ratio of $b_0^2/2\lambda =1$, so that the decay of the impulsive peak depends on the velocity $b_0$ as $I_{{\beta }^{*}}\left(t\right)\sim e^{-b_0^{2}t/2}$.

5. Quantum Mechanical Extension

There is an intertwined network of mutual analogies— between the theory of Brownian motion (and more generally stochastic dynamics) and quantum mechanics— that has led to fruitful physical insights and useful mathematical representations, borrowing ideas and techniques from one field to the other [49,50,51,52]. More specifically, there is a strong connection between stochastic processes that possess a finite propagation velocity and quantum mechanics whenever the relativistic constraint on the upper bound for velocities in the space-time propagation of physical phenomena is taken into account. This was already observed by Feynman [53] in his development of the path-integral formulation of quantum mechanics as regards the path of a relativistic free particle (see Figure 2.4 in [53]). However, Gaveau et al. [37] provided the derivation of the strict analogy between the $1+1$ Dirac equation (1 spatial dimension + 1 temporal dimension) and the Poisson–Kac processes in the presence of imaginary transition rates, i.e.,:

$$\begin{array}{ccc}\hfill \frac{\partial {\psi }_{+}(x,t)}{\partial t}& =& -c\frac{\partial {\psi }_{+}(x,t)}{\partial t}+i{\lambda }_0\left[{\psi }_{+}(x,t)-{\psi }_{-}(x,t)\right]\hfill \\ \hfill \frac{\partial {\psi }_{-}(x,t)}{\partial t}& =& c\frac{\partial {\psi }_{-}(x,t)}{\partial t}+i{\lambda }_0\left[{\psi }_{-}(x,t)-{\psi }_{+}(x,t)\right],\hfill \end{array}$$

(64)

where c is the velocity of light in vacuo, and $i=\sqrt{-1}$ and $({\psi }_{+},{\psi }_{-})$ are the two components of the vector-valued wave function. The spatial probability density is $|{\psi }_{+}{|}^2+{\left|{\psi }_{-}\right|}^2$, consistent with the Dirac theory [54]. As for the probabilistic theory of Poisson–Kac processes, Equation (64) converges towards the Schrödinger equation in the limit of $c,{\lambda }_0\to \infty$, keeping fixed the ratio $c^2/2{\lambda }_0$ (Kac limit). This model has been analyzed and generalized by several authors [55,56,57,58,59].

In this Section, we extend the stochastic model considered in the previous paragraphs to a quantum mechanical setting through the approach used by Gaveau et al. [37] discussed above, by considering a continuous velocity instead of a single velocity c, as in (64). The reason for this extension stems from the observation that the complete description of the dynamics of the stochastic process (43) requires the introduction of the addition variable $\beta$ parametrizing the partial density functions $p(x,t;\beta )$. In a quantum extension, this variable can be viewed as a “hidden variable” associated with a sub-quantum level of description of the system, from which classical quantum theory can be viewed as an emergent property. In other words, this quantum extension is aimed at providing an archetype for possible sub-quantum descriptions of reality, in the spirit of the work by David Bohm [38,39], which shows how classical quantum theory can be emergently derived from this “hidden” level. As we are primarily interested in highlighting these two aspects, the simple case of a free particle is considered in a one-dimensional spatial setting, in agreement with the stochastic model developed in the previous Sections.

The quantum mechanical counterpart of Equation (43) is given by the following equation:

$$\frac{\partial \psi (x,t;\beta )}{\partial t}=-b_0\beta \frac{\partial \psi (x,t;\beta )}{\partial x}+i{\lambda }_0\left[\psi (x,t;\beta )-\frac{1}{2}{\int }_{-1}^1\psi (x,t;{\beta }^{\prime })d{\beta }^{\prime }\right],$$

(65)

where $\psi (x,t;\beta )$ is the wavefunction at the subquantum level, parametrized with respect to the additional variable $\beta \in [-1,1]$. The constants $b_0$ and ${\lambda }_0$ are determined below.

Let $\Psi (x,t)$ be the classical (Schrödinger) wavefunction:

$$\Psi (x,t)={\int }_{-1}^1\psi (x,t;\beta )d\beta ,$$

(66)

then introduce the flux $J_{\psi }(x,t)=b_0{\int }_{-1}^1\beta \psi (x,t;\beta )d\beta$, such that:

$$\frac{\partial \Psi (x,t)}{\partial t}=-\frac{\partial J_{\psi }(x,t)}{\partial x}.$$

(67)

From Equation (65), it follows that $J_{\psi }(x,t)$ fulfills the following equation: Equation

$$\frac{\partial J_{\psi }(x,t)}{\partial t}=-b_0^2\frac{\partial }{\partial x}{\int }_{-1}^1{\beta }^2\psi (x,t;\beta )d\beta +i{\lambda }_0{J}_{\psi }(x,t).$$

(68)

Let us consider the Kac limit of this model, letting $b_0,{\lambda }_0\to \infty$ and keeping constant the nominal diffusivity $b_0^2/2{\lambda }_0=D_{\mathrm{nom}}$. From Equation (65), in the Kac limit, we have the following:

$$\psi (x,t;\beta )=\frac{1}{2}{\int }_{-1}^1\psi (x,t;\beta )d\beta =\frac{\Psi (x,t)}{2},$$

(69)

i.e., $\psi (x,t;\beta )$ in the Kac limit does not depend on $\beta$, while $J_{\psi }(x,t)$ attains the following expression:

$$J_{\psi }(x,t)=-i\frac{b_0^2}{2{\lambda }_0}2\frac{\partial }{\partial x}{\int }_{-1}^1{\beta }^2\psi (x,t;\beta )d\beta =-i\frac{2D_{\mathrm{nom}}}{3}\frac{\partial \Psi (x,t)}{\partial x},$$

(70)

where Equation (69) has been used. From Equation (70) inserted in Equation (67), it follows that if $D_{\mathrm{nom}}$ is such that $D_{\mathrm{nom}}=3D_{\hslash }/2$, where $D_{\hslash }=\hslash /2m$ is the quantum diffusivity for a massive particle of mass m, then $\Psi (x,t)$ satisfies the Schrödinger equation:

$$-i\frac{\partial \Psi (x,t)}{\partial t}=D_{\hslash }\frac{{\partial }^2\psi (x,t)}{\partial x^2}.$$

(71)

Moreover, it is easy to show that Equation (65) is a proper quantum mechanical equation, in the sense that a Born rule can be defined from it. In fact, elementary calculations provide the following:

$$\frac{\partial }{\partial t}{\int }_{-1}^1{\left|\psi (x,t;\beta )\right|}^{2}d\beta =-\frac{\partial }{\partial x}\left(b_0{\int }_{-1}^1\beta {\left|\psi (x,t;\beta )\right|}^{2}d\beta \right).$$

(72)

Therefore, ${\int }_{-1}^1{\left|\psi (x,t;\beta )\right|}^{2}d\beta$ is a conserved quantity corresponding to the quantum probability density function with respect to the position x. Observe that in the Kac limit, $\psi (x,t;\beta )\sim \Psi (x,t)$, and therefore ${\int }_{-1}^1{\left|\psi (x,t;\beta )\right|}^{2}d\beta \sim {|\Psi (x,t)|}^2$, which is in agreement with Schrödinger’s non-relativistic quantum theory.

For the sake of completeness, it is useful to observe that Equation (65) refers to a single free particle, but its extension to include the action of a potential or its generalization to a many-body problem is straightforward. If $U\left(x\right)$ is a potential acting on a massive particle, the associated quantum model generalizing Equation (65) is expressed by the following equation:

$$\begin{array}{ccc}\hfill \frac{\partial \psi (x,t;\beta )}{\partial t}& =& -b_0\beta \frac{\partial \psi (x,t;\beta )}{\partial x}-i\frac{U\left(x\right)}{\hslash }\psi (x,t;\beta )\hfill \\ & +& i{\lambda }_0\left[\psi (x,t;\beta )-\frac{1}{2}{\int }_{-1}^1\psi (x,t;{\beta }^{\prime })d{\beta }^{\prime }\right].\hfill \end{array}$$

(73)

Using the same approach outlined above, it is easy to see that Equation (65) converges in the Kac limit of the non-relativistic Schrödinger equation in the presence of the potential $U\left(x\right)$. In a similar way, the extension to a many-body problem is formally simple. To begin with, consider a single spatial dimension. Let $\mathbf{x}=(x_1,\dots ,x_n)$ be the coordinate vector of the particle system, where $x_k$ is the position of the k-th particle, $k=1,\cdots ,n$, and $\beta =({\beta }_1,\cdots ,{\beta }_n)$, $\lambda =({\lambda }_1,\cdots ,{\lambda }_n)$, so that $b_0{\beta }_k$ is the velocity of the k-th particle with ${\beta }_k\in [-1,1]$, and ${\lambda }_k$ is its transition rate. The extension of Equation (65) to an n-body problem on the real line in the presence of the interaction potential $U\left(\mathbf{x}\right)$ reads as follows:

$$\begin{array}{ccc}\hfill \frac{\partial \psi (\mathbf{x},t;\beta )}{\partial t}& =& -b_0\sum _{k=1}^n{\beta }_k\frac{\partial \psi (\mathbf{x},t;\beta )}{\partial x_k}-i\frac{U\left(\mathbf{x}\right)}{\hslash }\psi (\mathbf{x},t;\beta )\hfill \\ & +& i\sum _{k=1}^n{\lambda }_k\left[\psi (\mathbf{x},t;\beta )-\frac{1}{2^n}{\int }_{-1}^1\cdots {\int }_{-1}^1\psi (\mathbf{x},t;{\beta }^{\prime })d{\beta }^{\prime }\right],\hfill \end{array}$$

(74)

where $d{\beta }^{\prime }=d{\beta }_1^{\prime }·\cdots ·d{\beta }_n^{\prime }$. In higher dimensions, considering only translational motion, let us define $\mathbf{X}=({\mathbf{x}}_1,\cdots ,{\mathbf{x}}_n)$, where ${\mathbf{x}}_k=(x_{k,1},\cdots ,x_{k,d})$ is the position vector of the k-th particle in a d-dimensional space, and $\overline{\beta }=({\beta }_1,\cdots ,{\beta }_n)$, where the velocity of the k-th particle is $b_0{\beta }_k=b_0({\beta }_{k,1},\cdots ,{\beta }_{k,d})$. Thus, assuming isotropic transition rates, Equation (74) simply becomes:

$$\begin{array}{ccc}\hfill \frac{\partial \psi (\mathbf{X},t;\overline{\beta })}{\partial t}& =& -b_0\sum _{k=1}^n{\beta }_k·{\nabla }_k\psi (\mathbf{X},t;\overline{\beta })-i\frac{U\left(\mathbf{X}\right)}{\hslash }\psi (\mathbf{X},t;\overline{\beta })\hfill \\ & +& i\sum _{k=1}^{n}d{\lambda }_k\left[\psi (\mathbf{X},t;\overline{\beta })-\frac{1}{2^{nd}}{\int }_{-1}^1\cdots {\int }_{-1}^1\psi (\mathbf{X},t;{\overline{\beta }}^{\prime })d{\overline{\beta }}^{\prime }\right],\hfill \end{array}$$

(75)

where ${\nabla }_k=\left(\partial /\partial x_{k,1}\cdots \partial /\partial x_{k,d}\right)$ and $d{\overline{\beta }}^{\prime }={\prod }_{k=1}^n{\prod }_{h=1}^{d}d{\beta }_{k,h}^{\prime }$. The inclusion of n-body interactions expressed by the potential $U\left(\mathbf{X}\right)$ is perfectly admissible in a non-relativistic framework, while it would require care and attention from a relativistic perspective. But this discussion is manifestly outside the scope of the present work.

Spectral Properties and Dispersion Relations

Let us consider the spectral properties of the generalized Schrödinger operator,

$${\mathcal{L}}_q\left[\varphi (x,\beta )\right]=-b_0\beta \frac{\partial \psi (x,\beta )}{\partial x}+i{\lambda }_0\left[\psi (x,\beta )-\frac{1}{2}{\int }_{-1}^1\psi (x,{\beta }^{\prime })d{\beta }^{\prime }\right],$$

(76)

and let $\epsilon$ be an eigenvalue, ${\mathcal{L}}_q\left[\varphi (x,\beta )\right]=\epsilon \varphi (x,\beta )$. The eigenfunctions are of the form $\varphi (x,\beta )=e^{ikx}{\psi }_0\left(\beta \right)$. As before, by introducing the normalized quantities ${\epsilon }^{*}=\epsilon /{\lambda }_0$, $k^{*}=k{b}_0/{\lambda }_0$, it thus follows that that the “internal component” of the eigenfunction ${\varphi }_0\left(\beta \right)$ associated with the normalized eigenvalue ${\epsilon }^{*}$ is given by the following equation:

$${\varphi }_0\left(\beta \right)=\frac{A}{{\epsilon }^{*}-1+k^{*}\beta },$$

(77)

where A is a constant, leading to the following eigenvalue equation:

$$-2={\int }_{-1}^1\frac{d\beta }{{\epsilon }^{*}-1+k^{*}\beta },$$

(78)

which provides, upon quadraturae and algebraic manipulations, the expression for the eigenvalue:

$${\epsilon }^{*}=1-\frac{k^{*}(1+e^{-2k^{*}})}{1-e^{-2k^{*}}}.$$

(79)

In the quantum case, for any wave vector k, there exists an eigenvalue in the point spectrum of the operator ${\mathcal{L}}_q$, which is not the case for its probabilistic counterpart ${\mathcal{L}}_{\beta }$, as thoroughly analyzed in Section 4.1.

Let us consider the scaling properties of the eigenvalue spectrum, i.e., the dispersion relation connecting the particle energy $E=-\hslash \epsilon$ to the wave vector k. The Kac limit ($b_0,{\lambda }_0\to \infty$, $b_0^2/2{\lambda }_0=\mathrm{constant}$) corresponds to small values of $k^{*}$. Expanding the exponential displayed in Equation (79) in a Taylor series up to the leading order, in this case $e^{-2x}=1-2x+2x^2-4x^3/3$, and by developing the necessary algebra, one obtains the following from Equation (79):

$${\epsilon }^{*}=-\frac{{\left(k^{*}\right)}^2}{3}+O\left({\left(k^{*}\right)}^3\right).$$

(80)

Neglecting the $O\left({\left(k^{*}\right)}^3\right)$ term, it implies the following for the eigenvalue $\epsilon$:

$$\epsilon ={\lambda }_0{\epsilon }^{*}=-\frac{{\lambda }_0}{3}\frac{b_0^2{k}^2}{{\lambda }_0^2}=-\frac{2D_{\mathrm{nom}}}{3}{k}^2=-D_{\hslash }{k}^2,$$

(81)

consequently giving the energy as follows:

$$E=-\hslash D_{\hslash }{k}^2=\frac{{\hslash }^2{k}^2}{2m},$$

(82)

which is consistent with the classical quantum mechanical result.

Consider the quantity $-\epsilon /D_{\hslash }$, which in the classical quantum limit is equal to $k^2$, and set $\sigma ={\lambda }_0/D_{\hslash }$. Since $k^{*}=\sqrt{3/\sigma }k$ it follows that:

$$-\frac{\epsilon }{D_{\hslash }}=-\sigma {\left[\frac{(1-e^{2k^{*}})-k^{*}(1+e^{-2k^{*}})}{1-e^{-2k^{*}}}\right]}_{k^{*}=\sqrt{3/\sigma }k}.$$

(83)

For large values, $\left|k\right|$, $-\epsilon /D_{\hslash }=\sigma \left|k^{*}\right|$, and therefore $-\epsilon /D_{\hslash }$ admits the following crossover behavior:

$$-\frac{\epsilon }{D_{\hslash }}=\left\{\begin{array}{ccc}{k}^2,& & \left|k\right|<\sqrt{3\sigma }\\ \sqrt{3\sigma }\left|k\right|,& & \left|k\right|>\sqrt{3\sigma }\end{array}\right..$$

(84)

This phenomenon is depicted in Figure 13. The linear scaling with $\left|k\right|$ in the high wave vector region is the signature in the quantum model of the bounded propagation velocity characterizing the generalized Schrödinger operator ${\mathcal{L}}_q$ for finite values of $b_0$ and ${\lambda }_0$.

6. Concluding Remarks

In this article, we have analyzed the spectral (eigenvalue) properties of several classes of stochastic processes that possess a finite propagation velocity, showing the occurrence of a lower bound for the real part of the eigenvalue spectrum. As regards this spectral property, it is immaterial whether the long-term diffusive behavior is regular (linear Einsteinian scaling of the mean square displacement) or anomalous (superdiffusive). This result marks a fundamental and significant difference with respect to the parabolic diffusion model, for which short wavelength perturbations decay at a rate proportional to the squared norm of the wave vector.

This result has several general implications. (i) In the approach to the foundations of the thermodynamics of irreversible processes [60]: The results obtained confirm the analysis developed in [61], indicating that the Markov operators associated with the SEOs of stochastic processes characterized by a finite propagation velocity are invertible and form a group of transformations parametrized with respect to time t. Conversely, in the case of parabolic diffusion models, the weaker semi-group property for $t\ge 0$ holds. (ii) In classical field and transport theory: The results indicate that the hyperbolic field equations (e.g., for mass, heat, and momentum transfer), which can be developed starting from microscopic equations of motions expressed in the form of stochastic equations and characterized by a bounded velocity of propagation [47], may have completely different stability properties than their parabolic counterparts.

The occurrence of a lower bound for the real part of the eigenspectrum of operators describing the statistical evolution of stochastic processes possessing a finite propagation velocity also occurs for continuum models with a spectrum of particle velocities. In this article, we considered a prototypical example on the real line, in which the relaxational properties at high wave vectors are entirely associated with the essential part of the spectrum. The lower spectral bound is a consequence of the internal recombination dynamics of the stochastic states of the system and corresponds to the largest decay exponent (with reversed sign) of the overall probability density function $P(x,t)$, for an initial configuration with all the particles located at the same point and possessing the same velocity.

From a broader perspective, the analysis of SEOs outlined in this article can be extended to other relevant physical phenomenologies. This is the case for the linear response of stochastic systems that possess anomalous behavior, addressed in [62,63,64,65,66] by either using Continuous Time Random Walk scalings or approximate SEOs involving fractional derivative operators. The same problem can be framed within the age formalism of LWs introduced in Section 3, by keeping the same definitions and evolution in Equations (7)–(10), and modifying the boundary condition of Equation (11) in the following form:

$$p_{\pm }(x,t;0)=A_{\pm ,+}\left(t\right){\int }_0^{\infty }\lambda \left(\tau \right)p_{+}(x,t;\tau )d\tau +A_{\pm ,-}\left(t\right){\int }_0^{\infty }\lambda \left(\tau \right)p_{-}(x,t;\tau )d\tau ,$$

(85)

where the transition probabilities $A_{\pm ,\pm }\left(t\right)$, $A_{\pm ,\pm }\left(t\right)\ge 0$, $A_{+,\pm }\left(t\right)+A_{-,\pm }\left(t\right)=1$ account for the effect of the forcing field $f\left(t\right)$ and are defined according to [62,65] as follows:

$$A_{+,\pm }\left(t\right)=\frac{1+\mu f\left(t\right)}{2},A_{-,\pm }\left(t\right)=\frac{1-\mu f\left(t\right)}{2},$$

(86)

where $0<\mu <1$ and $\left|f\right(t\left)\right|\le 1$. In the analysis of this problem, the spectral results outlined in this manuscript are of limited use, as the simplest way to tackle this problem from physical grounds is to consider the partial moment hierarchy $m_{\pm }^{\left(n\right)}={\int }_{-\infty }^{\infty }{x}^n{p}_{\pm }(x,t;\tau )dx$, which will be addressed in forthcoming works.

Getting back to our probabilistic model with a spectrum of velocities, its quantum extension displays two main qualitative properties: (i) Contrary to the purely probabilistic case, for any $k\in \mathbb{R}$, an eigenvalue of the quantum operator ${\mathcal{L}}_q$ exists; (ii) The boundedness in the propagation velocity characterizing ${\mathcal{L}}_q$ implies, in the quantum case, a linear dispersive relation of the energy E with respect to $\left|k\right|$ for high values of k. This result is a direct consequence of the undulatory propagation, as from quantization $E=\hslash \omega$, and the $\omega$ is related to the group velocity $v_g$, by $d\omega /dk=v_g$, implying $E\sim k$. The quantum model analyzed in Section 4.1 is a simple example of a sub-quantum model possessing internal (“hidden”) degrees of freedom (in this case, the internal variable $\beta \in [-1,1]$) and providing the same emergent results as classical quantum theory.