A compounded random walk for space-fractional diffusion on finite domains

Many physical processes are known to exhibit anomalous diffusion where the mean squared displacement (MSD) scales nonlinearly with time. A standard approach to modeling this is through continuous-time random walks (CTRWs) Montroll and Weiss (1965); Metzler and Klafter (2000). The CTRW describes a random walk in which the time between jumps and the jump lengths are stochastic. In the appropriate continuum and diffusion limit, the CTRW master equation can model superdiffusion Metzler and Klafter (2000) where the MSD grows faster than linearly with time.

Superdiffusion is of interest across a wide range of physical systems. Examples include: random walks with spatial memory Metzler and Klafter (2000); Fedotov et al. (2018), reaction-diffusion equations for evanescent particles Abad et al. (2010), the motion of endogenous particles within amoeba Reverey et al. (2015), the trajectories of intracellular lipid bound vesicles Fedotov et al. (2018), diffusion in magnetic resonance imaging data Magin et al. (2013), electrical propagation in heterogeneous media such as biological tissues Bueno-Orovio et al. (2014), intrusion of contaminated water through sand Benson et al. (2000); Kelly and Meerschaert (2019), and in aquifers Yin et al. (2020); Sun et al. (2020), animal foraging strategies Humphries et al. (2010), laser driven dissipative plasmas Liu and Goree (2008), bacterial swarms Mukherjee et al. (2021) and disordered quantum systems Mildenberger et al. (2007).

Superdiffusion is traditionally generated by CTRWs where the jump length distributions have heavy tails and unbounded variance, namely Lévy flights Metzler and Klafter (2000). The heavy tailed jump lengths of Lévy flights are considered pathological to modeling applications as boundary conditions become problematic Dybiec et al. (2017). In physical systems, both boundary conditions and accounting for potential fields during a flight can be critical. For example, amoebae, intracellular vesicles, contaminated water molecules and foraging animals feel external potential fields generated by their heterogeneous environment during uni-directional flights. It would be unphysical if such random walks could leap over potential wells. The compounded random walk in this paper addresses the pathological nature of Lévy flights on bounded domains and in potential fields while maintaining the same defining characteristics on the unbounded domain.

For CTRWs with finite mean waiting times and Lévy flights for jump lengths, the governing equation limits to the space-fractional diffusion equation on the unbounded domain Compte (1996); Metzler and Klafter (2000)

where $\alpha\in(0,1]$ is the fractional exponent, $D_{\alpha}$ is the fractional diffusion coefficient and $\rho(x,t)$ is the probability density function (PDF) of finding the particle at position $x$ at time $t$. Here the fractional Laplacian $-(-\nabla^{2})^{\alpha}$ is defined via a Riesz derivative, which can be expressed in various ways Bayın (2016); Cai and Li (2019); Lischke et al. (2020), with the defining property from the Fourier transform

The particle motion governed by (1) is characterized by a pseudo superdiffusive MSD, $\lim\limits_{\delta\rightarrow 2\alpha}\langle|x|^{\delta}(t)\rangle\sim t^{\frac{1}{\alpha}}$, and $\delta<2\alpha$ Metzler and Klafter (2000). While the underlying CTRW stochastic process behind (1) is broadly accepted on infinite domains, there are numerous physical problems on finite domains for which there is no general consensus Lischke et al. (2020). Some critical questions arise when considering random walks in finite domains: How can a random walk take its next jump if the jump length exceeds the size of the domain? How can such a non-local random walk incorporate possible effects from spatially varying forces, including potential barriers? These questions are developing areas of research with pertinent examples being stopped Lévy flights Baeumer et al. (2018); Kelly et al. (2019) and subordinated Brownian motion Zoia et al. (2007); Garbaczewski and Stephanovich (2019); Lischke et al. (2020); Sokolov et al. (2001); Brockmann and Sokolov (2002). Currently, space-fractional models with the Riesz fractional derivative that incorporate zero flux boundary conditions are unphysical requiring ad-hoc exterior conditions Liu and Goree (2008); Garbaczewski and Stephanovich (2019); Brockmann and Sokolov (2002). This is further emphasized by the lack of consistent solutions for the space fractional diffusion equation defined using Riesz fractional derivatives Monroy and Raposo (2024).

Alternatively, the fractional order Laplacian can be defined in the spectral sense through the eigenequation Zoia et al. (2007); Lischke et al. (2020),

While the Riesz and spectral fractional Laplacians are equivalent on the infinite domain Zoia et al. (2007); Lischke et al. (2020), the spectral fractional Laplacian remains well defined on finite domains. In order to use this for modeling purposes it requires a physical framework such as an underlying random walk.

In this paper, we introduce a compounded random walk whose governing equation gives rise to the space-fractional diffusion equation (1) with the fractional Laplacian defined spectrally (3) in the diffusion limit. The compounding effect will model a process where a particle switches between a waiting process and another random walk on a faster time scale. The emergence of the spectral Laplacian from the underlying random walk provides a much needed physical interpretation of the space-fractional diffusion equation on bounded domains. Furthermore, this compounded random walk leads to a space-fractional Fokker-Planck equation on both bounded and unbounded domains with proper path sampling of position dependent forces. Our microscopic model is significant for two reasons: Firstly, it provides a method for modeling superdiffusive processes, which are ubiquitous, when a potential is applied to their domain. Secondly, it provides insight into the underlying dynamics for a system whose behavior is well described by a macroscopic spectral Fokker-Planck equation. We note some related work on this in Sokolov et al. (2001); Brockmann and Sokolov (2002), but the vast bulk of the literature on random walk models for space fractional Fokker-Planck equations has not considered the spectral Fokker-Planck equation. Our random walk interpretation is naturally suited for Monte Carlo simulation methods, which we verify by comparison to both analytical and numerical solutions. Next, we outline the compounded random walk and its governing equation through a CTRW and find spectral representations for the random walker’s probability mass function.

Consider a particle, performing a CTRW on a lattice, that waits an exponentially distributed amount of time at site $x_{i}$, with rate $\gamma$, before jumping to a new site $x_{j}$ governed by a jump distribution $\Lambda(x_{j},x_{i})$. The lattice is taken to be one-dimensional with spacing $\Delta x=(b-a)/(m-1)$ on a finite interval $[a,b]$ with $m$ lattice points. We denote $\rho(x_{i},t)$ as the probability of finding the particle at site $x_{i}=a+(i-1)\Delta x$ ($1\leq i\leq m$) at time $t$. Due to probability conservation, $\sum_{j}\Lambda(x_{j},x_{i})=1$. The master equation for the evolution of $\rho(x_{i},t)$ is Montroll and Scher (1973); Angstmann et al. (2015a)

Our aim is to find a jump distribution, $\Lambda(x_{i},x_{j})$, for a compounded random walk where a jump consists of a random number of Markovian steps. The steps, which make up a single jump, each locally sample the potential which governs its next step, while all occurring instantaneously so that no time passes over the duration of a single jump. The physical interpretation of this is that each jump is constructed from another CTRW with the same potential but with a negligible waiting time in comparison to compounded CTRW. Figure 1 shows an example of a compounded random walk. It is a necessity that after each of these steps the particle remains confined within the domain and thus have finite variance.

We begin by defining a nearest neighbor single step distribution. The particle can step one lattice site to the left or right except if the particle attempts to leave the domain ($1\leq i\leq m$), then it remains at its current site. This can be expressed as the single step transition probability,

where $\delta_{ji}$ is the Kronecker delta, $r(x_{i})$ is the probability to jump to the right at position $x_{i}$ and $l(x_{i})=1-r(x_{i})$. The probability of a step to the right can be related to a potential, $V(x_{i})$, and an inverse temperature scale, $\beta$, via a Boltzmann weight Henry et al. (2010)

For the compounded jump distribution $\Lambda(x_{i},x_{j})$, in (4), we allow the particle to instantaneously step a random number of times $K$ in a single jump such that

where $p_{K}(k)=\text{Prob}\{K=k\}$. The $k$-step distribution is defined recursively using (5) as

Substituting the jump distribution (7) into the master equation (4), we obtain

In order to express the right hand side of (9) spectrally, we first write the eigenequation for the single step distribution,

which defines $X_{n}$ as an eigenfunction of the single step distribution with eigenvalue $(1-\mu_{n})$. Then the $k$-step distribution will share the eigenfunctions from (10) such that,

Consequently, the eigenequation for the jump distribution (7) is

where

is the probability generating function (PGF), given that $p_{K}(0)=0$.

Now, we use the spectral property of the jump distribution (12) to find the solution of (9) by writing,

where for brevity we let $A_{K,n}=1-G_{K}(1-\mu_{n})$ and $c_{n}$ is a coefficient determined by the initial condition $\rho(x_{i},0)$. Differentiating (14) gives

Referring back to (4), it follows from the equivalence of the right hand sides of (9) and (15) that any linear spatial operator of the form

has an equivalent spectral representation of the form

where, as above, $A_{K,n}$ can be found from the generating function for the compounded random walk. The validity of (16) and (17) relies on the existence of $m$ distinct eigenvalues. This is ensured as the eigenvalue problem for the single step distribution (10) is equivalent to that for a real tridiagonal matrix with strictly positive values. The jump distribution (12) having $m$ distinct eigenvalues automatically follows as the PGF (13) is monotonically increasing.

Note that this spectral representation is general enough to account for any compounded nearest neighbor random walk in bounded domains and with position dependent potential fields. The master equation for the compounded random walk, (4) and (9), can be written in terms of the linear spatial operator (17) as

To recover the spectral Laplacian (3) from the discrete space random walk governed by (18), we choose a specific compounding distribution. As each of the steps have finite variance, to recover the superdiffusive properties of the random walk, we mandate that the distribution of the number of steps itself have infinite mean. For simplicity, we will use the Sibuya distribution defined as Sibuya (1979); Sibuya and Shimizu (1981)

where $0<\alpha\leq 1$, we find the PGF (13) gives

The compounded random walk with the Sibuya distribution represents a walk where the number of steps within a jump follows a modified Bernoulli process where the probability of stopping on the $k$^th step, conditional on not stopping before, is $\alpha/k$. While the Sibuya distribution was used to model power-law waiting times in discrete time random walks Angstmann et al. (2015b), here we use it to recover the power-law jump lengths as a sum of local steps, when the domain is unbounded.

Substituting the PGF of the Sibuya distribution (20) into the spectral representation of the linear spatial operator for the compounded random walk (17), we obtain

The linear operator (21) for the unbiased case where $r(x_{i})=1/2$ is a discrete space operator analogous to the spectral fractional Laplacian in (3). Next, we take the diffusion limit of (18) to obtain the fractional diffusion equation (1), where the linear spatial operator (21) reverts to the Riesz derivative (2) on unbounded domains. In addition, we show that the governing equation (18) gives rise to a space-fractional Fokker-Planck equation with a spectrally defined operator on both bounded and unbounded domains.

Extending $X_{n}(x_{i})$ to continuous arguments, $X_{n}(x)$, so that we can take a Taylor series, we substitute (5) and (6) into (10). The Taylor series of $X_{n}(x\pm\Delta x)$ about $\Delta x=0$ gives,

Using the Fokker-Planck operator,

we can rewrite (22) as

To obtain boundary conditions, we can substitute (5) and (6) into (10) at the boundary points $x_{1}=a$ and $x_{m}=b$ to obtain,

and

As $\Delta x\rightarrow 0$, we see that the eigenfunctions of the Fokker-Planck operator

must equate to the limit of (24) since $X_{n}\rightarrow\chi_{n}$ in the continuum limit. This equivalence between (24) and (27) gives the relationship

The eigenfunctions $\chi_{n}$ are restricted by the continuum limit of (25) and (26) which impose Robin boundary conditions

for $x=a$ and $x=b$.

Taking the diffusion limit of (21) and using $\mu_{n}=\Delta x^{2}\lambda_{n}/2$ from (28) as $\Delta x\rightarrow 0$, we obtain

where the mean waiting time is defined by $\tau=1/\gamma$ and

The master equation (18) defined via the linear spatial operator in continuous space (30) becomes,

From (32), we arrive at a space-fractional Fokker-Planck equation,

where $(-\mathcal{L})^{\alpha}$ is the spectral fractional Fokker-Planck operator. From the eigenequation (27), the spectral fractional Fokker-Planck operator is defined by

The significance of (33) lies in its connection with an underlying compounded random walk whose $K$ random steps in each jump is Sibuya distributed (19). Note that the fractional Fokker-Planck operator in (33) is fundamentally different to previous formulations that arose from Lévy flights or the Langevin equation on the unbounded domain Fogedby (1994, 1998); Jespersen et al. (1999); Metzler and Klafter (2000). To demonstrate this, we consider the Fourier transform of (33), defined on the unbounded domain

where we have set a constant drift, $2\beta dV/dx=v$ in the Fokker-Planck operator (23). In previous models, the fractional Fokker-Planck equation for constant drift $v$ was defined in the Fourier space by Jespersen et al. (1999); Benson et al. (2000); Meerschaert and Tadjeran (2004); Shen et al. (2008)

Equation (36) is also known as the fractional advection-dispersion equation Benson et al. (2000); Meerschaert and Tadjeran (2004); Shen et al. (2008); Wang and Barkai (2020). Comparing (35) and (36), it is clear that when $v=0$, the two formulations are equivalent and we obtain the space-fractional diffusion equation (1). This is expected as the Riesz fractional Laplacian is equivalent to the spectral fractional Laplacian for unbounded domains Lischke et al. (2020). However, a clear difference between (35) and (36) is that, in (35) the advection is felt across the entire path of the underlying random walk and thus encapsulated within the term $(|k|^{2}+ikv)^{\alpha}$, as opposed to only accounting for it after landing on a site as evident in the separate terms of $D_{\alpha}|k|^{2\alpha}+ikv$. Our fractional Fokker-Planck equation with the spectral fractional Fokker-Planck operator, (35), is mathematically consistent with Sokolov et al. (2001); Brockmann and Sokolov (2002). Figure 2 illustrates the difference between the two governing equations on the unbounded domain at four different times. Unlike the solution to the fractional advection-dispersion equation, (36), the solution to the spectral fractional Fokker-Planck equation, (35), does not remain symmetric about the peak. This is due to the compounding of the bias effect in the solution to (35).

For bounded domains, (33) is identical to the fractional diffusion equation defined via the spectral fractional Laplacian with Neumann boundary conditions when $dV/dx=0$ in (29). For a spectral fractional Laplacian formulation with Dirichlet boundary conditions, one needs to reformulate (5) and take the diffusion limit again, in which case we obtain (33) with eigenfunctions satisfying boundary conditions $\chi_{n}(a)=\chi_{n}(b)=0$.

Physically, the Dirichlet boundary conditions model the probability distribution for a particle at location $x$ at time $t$ that has never reached the absorbing boundary. This immediately provides a clear use for finding first-passage time distributions for CTRWs making the governing equation well motivated.

The underlying compounded discrete-space random walk for (33) with zero flux boundary conditions gives a Monte Carlo scheme that can be performed as follows:

1. 1.

   Partition the interval $[a,b]$ into $m-1$ segments of length $\Delta x=(b-a)/(m-1)$ such that there are $m$ nodes approximating position with the first node at $a$.

2. 2.

   Initialize a particle at some initial position $X_{0}=x_{i}$ and set initial time as $t=0$.

3. 3.

   Draw a random time $\Delta t\sim\text{Exp}(\gamma)$.

4. 4.

   Draw a Sibuya distributed random variable $K$ Hofert (2011).

5. 5.

   Given $X=x_{j}$, update $X\leftarrow X+\xi\Delta x$ with the following rules:

   1. a)

      If $1<j<m$, draw $\xi$ with probabilities $\text{Prob}(\xi=1)=r(X)$ and $\text{Prob}(\xi=-1)=1-r(X)$;

   2. b)

      If $j=1$, then $\text{Prob}(\xi=0)=1-r(X)$ and $\text{Prob}(\xi=1)=r(X)$;

   3. c)

      If $j=m$, then $\text{Prob}(\xi=0)=r(X)$ and $\text{Prob}(\xi=-1)=1-r(X)$.

6. 6.

   Repeat step (5) an additional $K-1$ times.

7. 7.

   Update $t\leftarrow t+\Delta t$.

8. 8.

   Repeat steps (3)-(6) until $t\geq T$, where $T$ is the simulation end time.

The time discretization of (9) gives a numerical scheme for the solution of (14).

Figure 3 shows the Monte Carlo simulation, numerical approximation and the exact spectral solution to (33) with $dV/dx=0$, on the domain $x\in[-1,1]$ with Neumann boundary conditions and initial condition $\rho(x,0)=\delta(x)$. The analytical solution is

Now, we extend our results to the biased CTRW with a constant forcing term $\beta V^{\prime}(x)=-2$. The analytical solution for the governing equation (33), given initial condition $\rho(x,0)=\delta(x)$, is

where $\omega_{n}=n\pi/2$

and

Figure 4 shows the correspondence between the Monte Carlo simulation, numerical approximation and the truncated spectral solution (38) for the fractional spectral Fokker-Planck operator with constant bias in (33). As the solution evolves over time, it is evident that it approaches the Boltzmann distribution, a property not present in other random walk formulations Sokolov et al. (2001).

The spectral fractional Fokker-Planck equation may be used to model the survival probability in the first-passage time problem. The solution to (33) in the case with no bias, absorbing boundaries on the domain $[-1,1]$ and initial condition given by the Dirac delta function $\rho(x,0)=\delta(x)$ is written as

This solution represents the probability the particle is located at $x$ at time $t$ and has not reached the boundary. Alternatively, if the particles were only absorbed when they had finished a compounded jump outside the domain then the process becomes analogous to killed Lévy flights on the bounded domain and thus its governing equation gives a finite difference approximation to the space-fractional diffusion equation defined with the Riesz derivative and zero exterior conditions $\rho(x,t)=0$ for $x\notin[-1,1]$ Zoia et al. (2007). Figure 5 shows the difference in behavior over time, where it can be seen that there is a significant difference between the case when particles can be absorbed during their compounded jumps (left) and the case where particles can only be absorbed at the end of the compounded jumps (right). This effectively demonstrates that the governing equation with the spectral Laplacian obtained from the compounded CTRW on the bounded domain models the first passage time of a process different to that of CTRW models incorporating Lévy flights.

Our results connect an underlying compounded random walk with the spectral fractional Laplacian in the diffusion limit. This provides a novel and physically valid model to derive the space-fractional diffusion equation on bounded domains. Furthermore, this connection naturally generates a spectral fractional Fokker-Planck equation (33) that is well defined on bounded domains and samples space-dependent forces at every point along a jump. The derivation of the fractional Fokker-Planck equation enables a description of random walkers that experience space-dependent forces over the entire path of a single heavy-tailed jump, improving previous formulations Fogedby (1994, 1998); Jespersen et al. (1999). In doing so, we open an avenue for models with heavy-tailed jumps where the random walk is sampling the environment throughout the jump and is confined to a finite region in space. Moreover, through the PGF, (13), new generalized operators based on experimentally observed microscopic mechanisms can be formulated using (12), (16) and (17) for models of transport, such as intracellular transport Reverey et al. (2015); Fedotov et al. (2018), foraging animals Humphries et al. (2010), and disorderd quantum systems Mildenberger et al. (2007). The governing equations derived through this approach will have convenient spectral solutions and a natural Monte Carlo simulation scheme.