On subdiffusive continuous time random walks with stochastic resetting

TL;DR

This paper investigates two models of subdiffusive processes with stochastic resetting, deriving exact statistical properties and exploring their potential applications in biological data analysis.

Contribution

It introduces and analyzes two distinct models of subdiffusion with resetting, highlighting their non-Markovian differences and providing exact solutions for key statistical measures.

Findings

Derived exact moments and stationary distributions for both models

Identified differences due to non-Markovian nature of subdiffusion

Extended models to include external forces like constant drift

Abstract

We analyze two models of subdiffusion with stochastic resetting. Each of them consists of two parts: subdiffusion based on the continuous-time random walk (CTRW) scheme and independent resetting events generated uniformly in time according to the Poisson point process. In the first model the whole process is reset to the initial state, whereas in the second model only the position is subject to resets. The distinction between these two models arises from the non-Markovian character of the subdiffusive process. We derive exact expressions for the two lowest moments of the full propagator, stationary distributions, and first hitting times statistics. We also show, with an example of a constant drift, how these models can be generalized to include external forces. Possible applications to data analysis and modeling of biological systems are also discussed.

Equations

$$W_0(x,s|x_0)=\frac{s^{\frac{\mu }{2}-1}}{2\sqrt{D}}{e}^{-\sqrt{\frac{s^{\mu }}{D}}|x-x_0|},$$

$${\partial }_t{W}_0(x,t|x_0)={}_{0}D{}_t^{1-\mu }D{\partial }_{xx}^2{W}_0(x,t|x_0),$$

$$⟨X^2(t)⟩=\frac{2D}{\Gamma (\mu +1)}{t}^{\mu }.$$

$$W_0(k,s|x_0)=\frac{\Psi (s)e^{ik{x}_0}}{1-\psi (s)g(k)}.$$

$${\partial }_t{G}_0(x,t|x_0)={}_{0}D{}_t^{1-\mu }D{\partial }_{xx}^2{G}_0(x,t|x_0),$$

$$G_0(x,s|x_0))=\frac{s^{\mu /2-1}}{2\sqrt{D}}\left(e^{-\sqrt{s^{\mu }/D}|x-x_0|}-e^{-\sqrt{s^{\mu }/D}|x+x_0|}\right).$$

$$S_0(s)=\underset{0}{\overset{\infty }{\int }}\mathrm{d}x{G}_0(x,s|x_0)=\frac{1}{s}\left(1-e^{-\sqrt{s^{\mu }/D}|x_0|}\right),$$

$${\rho }_0(s)=1-s{S}_0(s)=\exp \left(-\sqrt{s^{\mu }/D}|x_0|\right)$$

$$W(x,t|x_0)=e^{-rt}{W}_0(x,t|x_0)+\underset{0}{\overset{t}{\int }}d\tau r{e}^{-r\tau }{W}_0(x,\tau |0),$$

$$W(x,s|x_0)=W_0(x,r+s|x_0)+\frac{r}{s}{W}_0(x,r+s|0).$$

$$W(k,s|x_0)=\frac{\Psi (r+s)\left(e^{ik{x}_0}+r/s\right)}{1-\psi (r+s)g(k)}.$$

$$\left(D(s+r{)}^{1-\mu }{k}^2+s+r\right)W(k,s|x_0)=e^{ik{x}_0}+\frac{r}{s}.$$

$${\partial }_{t}W(x,t|x_0)=D{e}^{-rt}{}_{0}D{}_t^{1-\mu }{e}^{rt}{\partial }_{xx}^{2}W(x,t|x_0)-rW(x,t|x_0)+r\delta (x).$$

$$f_s(x)=\underset{t\to \infty }{\lim }W(x,t|x_0)=\underset{s\to 0}{\lim }sW(x,s|x_0)=\frac{1}{2}\sqrt{\frac{r^{\mu }}{D}}{e}^{-\sqrt{\frac{r^{\mu }}{D}}|x|}.$$

$$⟨X^2(s)⟩=\frac{2D}{s(s+r{)}^{\mu }}.$$

$$\text{mbox}VarX(t)=⟨X^2(t)⟩=\frac{2D}{\Gamma (\mu )}\underset{0}{\overset{t}{\int }}d\tau {\tau }^{\mu -1}{e}^{-r\tau }=\frac{2D}{r^{\mu }\Gamma (\mu )}\gamma (\mu ,rt),$$

$$⟨X^2(t)⟩\approx \frac{2D}{r^{\mu }}\left(1-\frac{e^{-rt}}{\Gamma (\mu )(rt{)}^{1-\mu }}\right),$$

$$\begin{array}{rl}{T}_r& =I[T_0\le R]T_0+I[T_0>R](R+T_r^{\prime })\\ & =\min (T_0,R)+I[T_0>R]T_r^{\prime },\end{array}$$

$$⟨T_r⟩=\frac{⟨\min (T_0,R)⟩}{⟨I[T_0\le R]⟩},$$

$$⟨T_r⟩=\frac{S_0(r)}{1-r{S}_0(r)},$$

$$⟨T_r⟩=\frac{1}{r}\left(e^{\sqrt{r^{\mu }/D}|x_0|}-1\right).$$

$$\underset{r=\frac{1}{\mu }\to \infty }{\lim }⟨T_r⟩=0,$$

$${\partial }_t⟨X^2(t)⟩\propto t^{\mu -1}.$$

$$r^{\ast }={\left(\frac{z_{\mu }^{2}D}{x^2}\right)}^{\frac{1}{\mu }},$$

$$\frac{\mu }{2}{z}_{\mu }=1-e^{-z_{\mu }}.$$

$$z_{\mu }\approx \frac{2}{\mu },$$

$$r^{\ast }\propto {\left(\frac{2}{\mu }\right)}^{\frac{2}{\mu }}.$$

$$⟨\min (T_0,R{)}^2⟩=-2{\partial }_r{S}_0(r)$$

$$⟨I[T_0>R]\min (T_0,R)⟩=-r{\partial }_r{S}_0(r),$$

$$\mathrm{Var}{T}_r=⟨T_r^2⟩-⟨T_r{⟩}^2=-\frac{2{\partial }_r{S}_0(r)+S_0(r{)}^2}{1-r{S}_0(r)}.$$

Full text

I Introduction

Recently a simple model of diffusion with stochastic resetting has been proposed majumdar2011resetting . In its basic form it is a Brownian motion interrupted by reset events, which instantaneously bring the process back to the initial position. Resets happen randomly in time, according to the Poisson point process with intensity $r$. In contrast to standard, symmetric diffusion, diffusion with resetting leads to a finite mean first arrival time (MFAT) to any fixed position majumdar2011resetting . Moreover, it has been shown that resetting may be beneficial for search in a variety of modified scenarios, e.g. when the initial random walk process is superdiffusive kusmierz2014first ; kusmierz2015optimal , when the space is higher-dimensional evans2014diffusion and/or bounded christou2015diffusion , and when the resetting intensity is inhomogeneous in time Reuveni2016 ; pal2016first ; eule2016non ; bhat2016stochastic ; pal2016diffusion ; kusmierz2018robust ; maso2019stochastic ; maso2019transport . The combination of a random walk process with stochastic resetting events leads to many nontrivial features: A non-local character of resets results in a current-carrying non-equilibrium steady states (NESS), i.e. stationary distributions in which detailed balance does not hold majumdar2011resetting ; evans2013optimal ; pal2015diffusion ; eule2016non ; montero . When analyzed as a search process involving a mixture of local steps with resets, the mean first passage time to find a target has been computed exactly and has been found to have a minimum at an optimal resetting rate majumdar2011resetting , indicating that the search process is more efficient in the presence of resetting. This non-equilibrium search is more efficient than a corresponding (leading to the same shape of the stationary distribution) conservative-force-induced equilibrium dynamics evans2013optimal and in some cases more efficient than any conservative-force-induced equilibrium dynamics kusmierz2017optimal . When instead of the Brownian motion the search is performed by long range Lévy moves combined with resets, the motion generates a bifurcation in the optimal search strategy (i.e. bifurcation in the shortest mean first-time to reach a target). The bifurcation is discontinuous kusmierz2015optimal or continuous kusmierz2014first , depending on whether the process is performed in discrete or continuous time. Stochastic resets can modify the splitting probabilities leading to a higher success rate in Bernoulli trials belan2018restart . A one-dimensional version of the walk where the walker, at each time step, resets to the maximum of the already visited positions with a certain probability leads to a dynamical transition in the temporal relaxation majumdar2015dynamical . The stochastic resetting model is related to reinforced random walks davis1990reinforced , a version of which has recently been proposed as a model of large-scale movement of free-ranging capuchin monkey and shown to lead to super-slow (logarithmic) diffusion boyer2014random . It can also lead to a localization transition similar to Anderson localization if supplemented by a simple model of non-homogeneous learning falcon2017localization .

Here we explore another generalization of diffusion with resetting, wherein the diffusive process between resets is substituted with a process generated by the continuous-time random walk (CTRW) scheme montroll1965random ; barkai2000continuous ; meerschaert2004limit ; metzler2000random . In contrast to recent works that combine stochastic resets with the CTRW with exponentially distributed waiting times between jumps montero ; campos2015phase ; Montero2017 , we assume that the waiting times are distributed according to a power law $p(t)\propto t^{-\mu-1}$, where $t\gg 1$ and $\mu\in(0,1)$. Due to the infinite mean waiting time the primary process (i.e. without resetting) is subdiffusive with the mean square displacement (MSD) growing sublinearly with time, $\langle X^{2}(t)\rangle\propto t^{\mu}$, and the Brownian diffusion is recovered in the limit $\mu\to 1$. A subdiffusive CTRW process is non-Markovian scalas2000fractional and shows weak egodicity breaking bel2005weak ; lubelski2008nonergodicity , i.e. time and ensemble averages of physical observables do not coincide.

In addition to the CTRW process there are several distinct mathematical models of subdiffusion sokolov2015 , differing significantly in the details of dynamics. Two most popular are fractional Brownian motion (FBM) mandelbrot1968fractional ; taqqu1975weak ; biagini2008stochastic ; mishura2008stochastic ; deng2009ergodic and diffusion on fractals d1983random ; o1985diffusion ; ben2000diffusion . Each of these models corresponds to a different physical mechanism leading to the subdiffusive behavior. While all of them share the sublinear growth of the variance, they exhibit many dissimilarities in other features. Importantly, in contrast to the subdiffusive CTRW, the other two processes are ergodic (but see deng2009ergodic ). Recently, statistical methods of discriminating between FBM and CTRW magdziarz2009fractional , diffusion on fractals and CTRW condamin2008probing , and diffusion on fractals and FBM meroz2013test have been proposed.

Subdiffusion has been observed in a number of distinct physical phenomena scher1975a ; kopelman1984reaction ; scher2002dynamical ; weiss2007anomalous . Many experiments have shown that the transport in living cells is subdiffusive schwille1999molecular ; wachsmuth2000anomalous ; tolic2004anomalous , but the underlying mechanism is not yet fully understood. In this context, some recent works report non-ergodic CTRW-like behavior jeon2011vivo ; weigel2011ergodic ; jeon2012anomalous ; tabei2013intracellular , whereas others are in favor of ergodic FBM or diffusion on fractals szymanski2009elucidating ; weber2010bacterial ; kepten2011ergodicity ; barkai2012single . In view of those findings, we study how the statistical properties of a subdiffusive CTRW process are altered under random relocations. The overall process is a combination of a CTRW and independent resets. In subsequent Sections we define two models of the resetting mechanism when either relocation event erases memory about the past motion (Model I) or, otherwise resets do not interfere with waiting for the next step of the random motion (Model II). As we will show, non-trivial features of such processes may not only provide a descriptive model of natural phenomena, such as dynamics of molecular motors in the crowded cytoplasm, but also may be helpful in designing resampling-based statistical inference methods aiming at detecting trapping and associated ergodicity breaking.

Some of our results can be considered to be special cases of previous more general considerations pal2016first ; shkilev2017continuous ; maso2019transport . In contrast to these works, we focus almost exclusively on the CTRW subdiffusion with exponential resetting and analyze the possible outcomes in depth, uncovering intriguing effects that, to best of our knowledge, have not yet been reported. Power-law waiting times have also been recently explored in evans2018effects ; maso2019stochastic in the context of refractory periods/residence times that follow resets. Such trapping at the resetting position is related to, but distinct from the space-homogeneous trapping which happens in our Model I. Interestingly, our Model II also has a net effect of trapping at the resetting position. Notice, however, that the underlying mechanism is different from models in evans2018effects ; maso2019stochastic and is driven by the power-law distributed waiting times between steps.

Throughout the paper we use the following notation. The propagator of a process of interest is written as ${W(x,t|x_{0})\equiv W(x,t|x_{0},0)}$, with the initial condition $W{(x,0|x_{0})=\delta(x-x_{0})}$. Partial derivatives with respect to time and space are written as $\partial_{t}$ and $\partial_{x}$, respectively. The parameter $r$, called resetting intensity, gives the expected number of reset events in unit time. We denote the space Fourier transform and the time Laplace transform implicitly, by changing the argument(s) of the transformed function ($x\to k$ and $t\to s$, respectively).

II Subdiffusion and Fractional Fokker-Planck equation

Here we list basic results related to CTRW subdiffusion that we will use throughout the paper. For a complete description, see e.g. metzler2000random . In our derivations of quantities related to the process with resetting we will often employ the function

[TABLE]

which is the (Laplace transformed) solution of a time-fractional Fokker-Planck (FFPE) equation

[TABLE]

where ${}_{0}D{{}_{t}^{1-\mu}}$ denotes the Riemann-Liouville fractional derivative operator metzler1999anomalous ; metzler2000random ; metzler2000subdiffusive with $0<\mu<1$. Function $W_{0}(x,t|x_{0})$ is the propagator of the process without resetting, i.e. probability density function (PDF) of the CTRW subdiffusion starting from $x_{0}$ at time $t_{0}=0$. Multiplying (1) by $(x-x_{0})^{2}$ and integrating over $x$ we obtain the variance (in this case equivalent to the MSD) evolution in the Laplace space $\langle X^{2}(s)\rangle=2D/s^{1+\mu}$ which translates to

[TABLE]

The FFPE (2) can be derived as follows. The movement of a particle is generated via jumps generated from a PDF $g(x)$. The time between two consecutive jumps is generated from a non-negative random variable, called waiting times, that is described by a PDF $\psi(t)$ or, equivalently, by the probability that the time before the next jump happened is longer than $t$, $\Psi(t)=\int_{t}^{\infty}d\tau\psi(\tau)$. Between subsequent jumps the particle is immobile and, importantly, jumps and waiting times are independent. It is well known that the propagator of such process is given by the Montroll-Weiss formula montroll1965random ; metzler2000random ; zaburdaev2015levy

[TABLE]

In the case of the subdiffusive CTRW, we assume that ${g(k)\approx 1-D_{0}k^{2}}$ and ${\psi(s)\approx 1-A_{\mu}s^{\mu}}$ and perform the limit ${A_{\mu}\to 0}$, ${D_{0}\to 0}$ with ${D_{0}/A_{\mu}\to D}$, which leads to (1) and (2). Most results obtained in this paper will have two forms: a general CTRW form generalizing (4) and a subdiffusive form that follows from taking the described limit. In some cases this will allow us to write a modified fractional Fokker-Planck equation. Nevertheless, our derivations are based on the renewal approach that can be used in the very general setting of the CTRW.

It is easy to calculate the Laplace transform of the survival probability $S_{0}(t)$ of the subdiffusive search process without resetting in the presence of a single target located at $x=0$. The corresponding unnormalized density $G_{0}(x,t|x_{0})$ is the solution of the equation

[TABLE]

which has the same form as (2) except for the boundary conditions. The target is introduced by imposing the absorption condition ${G_{0}(x=0,t|x_{0})=0}$. We can construct such a solution from the solution (1) of (2) via the method of images, leading to

[TABLE]

The corresponding survival probability reads

[TABLE]

which shows that, despite the trapping, in the lomg time limit $\lim_{t\to\infty}S_{0}(t)=0$, i.e. the target is almost surely (with probability $1$) reached. The first arrival times (FAT) probability density function

[TABLE]

behaves like $\propto t^{-(1+\mu/2)}$ for large $t$, which generalizes the well-known Sparre Andersen scaling $\propto t^{-3/2}$, which holds for symmetric Markovian processes with stationary increments andersen1954fluctuations ; andersen1955fluctuations ; dybiec2009anomalous .

III First model: resetting of the process

Resets in the first model bring the position to the origin and delete all memory about the past. In the CTRW scheme waiting times for the diffusive jump and for the reset are drawn independently. Whenever the latter is shorter, the position is reset to [math] and both waiting times are generated again. Thus, the aging effect of subdiffusion is also subject to reset. In the following we will refer to such events as hard resets. Although we assume that resetting brings the process back to the initial position, it is convenient to perform calculations for any initial position $x_{0}$, with the resetting position kept at the origin. At the end we set $x_{0}=0$.

Given these assumptions the propagator of the process subject to random resets can be written in terms of the propagator of the process without resetting via renewal equation (gupta2014fluctuating ; kusmierz2015optimal )

[TABLE]

which in the Laplace space reads

[TABLE]

We can therefore use the known propagator of the CTRW process (4) to construct the propagator of the process with hard resets, which reads

[TABLE]

Expression (11) is very general and can be used to study any CTRW process with hard resets. In fact, it would be straightforward to generalize it further to include spatio-temporal coupling between the jump distribution and waiting times, as in the case of Lévy walks zaburdaev2015levy . As announced, here we focus on the special case of subdiffusive CTRW.

We note in passing that, for any $g(x)$ that is continuous at $x=0$, the propagator (11) has a singular part that behaves like $\Psi(r+s)(e^{ikx_{0}}+r/s)$. However, if we assume that ${g(k)\approx 1-D_{0}k^{2}}$ and ${\psi(s)\approx 1-A_{\mu}s^{\mu}}$ with ${A_{\mu}\to 0}$, ${D_{0}\to 0}$, and ${D_{0}/A_{\mu}\to D}$, the singular part disappears. As we will show later, this is one of the features that discriminates hard from soft (Model II) resets.

Taking such limit is equivalent to plugging the subdiffusive $W_{0}$ given by (1) into (10).

The corresponding propagator with hard resets in the Fourier-Laplace space $W(k,s|x_{0})$ solves the following algebraic equation

[TABLE]

Eq. (12) corresponds to the following integro-differential equation in $(x,t)$ space

[TABLE]

The operator $e^{-rt}{}_{0}D{{}_{t}^{1-\mu}}e^{rt}$ corresponds to a version of truncated Riemann-Liouville fractional operator meerschaert2008tempered ; gajda2010fractional ; alrawashdeh2017applications . The same operator was previously proposed in the context of reaction-subdiffusion schmidt2007mesoscopic and, interestingly, it also appears in the context of subdiffusion with a constant external force that affects waiting times fedotov2015subdiffusion . In our case, it emerges because the resets induce two effects: the position of the particle is set to [math] and the waiting times are generated anew. The latter phenomenon introduces an effective truncation on the waiting times, which is realized in the Fractional Fokker-Planck equation through the described modification of the Riemann-Liouville fractional operator. From Eq. (12) we obtain the stationary distribution

[TABLE]

Recently the same formula has been independently derived in maso2019transport . As expected, for $\mu=1$ we recover formula derived for diffusion with stochastic resetting majumdar2011resetting . For any $\mu$ we obtain the Laplace distribution. The difference lies in the dependence of its mean value on the resetting intensity, which could also be deduced from dimensional analysis.

In order to examine the transient behavior of the process we calculate first and second moments. First, we focus on the general case of a process starting from some position $x_{0}$. By inspecting Eq. (12) it is easy to check that $\langle X(t)\rangle=x_{0}e^{-rt}$, i.e. relaxation of the expected position is exponential and is driven purely by the resetting events. Relaxation of the variance incorporates two distinct processes. The first of them is related to relaxation of the position, and is described by the expression $(e^{-rt}-e^{-2rt})x_{0}^{2}$. As announced before, we focus on the second relaxation process, which is strictly related to the width change, and thus we set $x_{0}=0$, i.e. we assume that the process starts from the resetting position. From Eq. (12) we obtain the MSD in the Laplace space

[TABLE]

Inverting the Laplace transform leads to the integral expression

[TABLE]

where $\gamma(s,x)$ stands for the lower incomplete gamma function. Asymptotically

[TABLE]

as $t\to\infty$, i.e. the stationary state is approached exponentially fast.

In the following we calculate the first arrival time (FAT) statistics. We follow the general strategy introduced in pal2016first . Due to the renewal structure of the resetting, we can link the FAT distribution of the process with resetting with $S_{0}$. We introduce a random variables $R$ (time to the next reset), $T_{0}$ (FAT of the process without resetting), and $T_{r}$ (FAT of the process with resetting). The resets can be introduced as follows: first draw $T_{0}$ and $R$. If $T_{0}\leq R$ then the reset was supposed to happen after the target was found, thus $T_{r}=T_{0}$. If, on the other hand, $T_{0}>R$, then the reset came before the target was supposed to be found, thus the search has to start anew. In that case $T_{r}=R+T_{r}^{\prime}$, where $T_{r}^{\prime}$ is a new realization of the FAT of the process with resetting, which in turn depends on the independent realizations of the resetting times and FATs of the process without resetting. We can summarize this procedure with the equation

[TABLE]

where the indicator function ${I[a>b]=1}$ if ${a>b}$ and ${I[a>b]=0}$ if ${a\leq b}$. Since $T_{r}^{\prime}$ is independent from $T_{0}$ and $R$, we can easily average both sides of (18) leading to the known general formula for the MFAT of the process subject to resetting pal2016first

[TABLE]

where averaging is performed over independent $R$ and $T_{0}$. In the case of exponential (i.e. constant rate) resetting ${\langle\min(R,T_{0})\rangle=S_{0}(r)}$ and ${\langle I[T_{0}\leq R]\rangle=1-rS_{0}(r)}$, thus

[TABLE]

where $S_{0}(r)\equiv S_{0}(s=r)$. Combining (7) and (20) we arrive at

[TABLE]

Equation (21) is consistent with a more general formula derived with different methods in shkilev2017continuous . In the limit of $\mu\to 1$ it simplifies to the well-known formula for the MFAT of diffusion with stochastic resetting majumdar2011resetting . In contrast to superdiffusive Lévy flights with resetting kusmierz2014first ; kusmierz2015optimal , the CTRW subdiffusion with resetting preserves the exponential form of of the dependence of the MFAT on the distance to the target. Interestingly, for a fixed $|x_{0}|>0$, there is no local minimum in parameter space $(r,\mu)$. The infimum of the MFAT is [math],

[TABLE]

because subdiffusion is slow at large time-scales, but extremely fast at short time-scales. This can be seen from the following formula, which holds for pure subdiffusion without resetting ($r=0$)

[TABLE]

Note, however, that at short time-scales this formula cannot describe physical phenomena for which the speed is bounded.

For a fixed $\mu$ there is exactly one minimum of the MFAT as a function of $r$. The optimal resetting intensity is given by the expression

[TABLE]

where $z_{\mu}$ is the solution of the transcendental equation

[TABLE]

When $z_{\mu}$ is large, the exponential element on the RHS of Eq. (25) is negligible, thus we can write an approximate solution

[TABLE]

which is valid for $\mu\ll 1$. This leads to the following dependence of the optimal resetting intensity on $\mu$

[TABLE]

How is $T_{r}$ distributed around its mean value? We square both sides of (18) and calculate averages

[TABLE]

and

[TABLE]

arriving at

[TABLE]

Combining (7) and (30) we obtain

[TABLE]

where $y=r^{\mu/2}|x_{0}|/\sqrt{D}$. Let us analyze the squared coefficient of variation.

[TABLE]

At $r=r^{*}$ (and thus $y=z_{\mu}$), i.e. when the intensity of resetting is optimal, the coefficient of variation is equal to $1$ (see Fig. 2), in line with the universal result obtained recently Reuveni2016 . The resetting intensity at which CV admits its minimum is given by

[TABLE]

where $y_{\mu}$ is the solution of the transcendental equation

[TABLE]

Note that $r^{**}>r^{*}$ for any $\mu\in(0,1]$ (cf. Fig. 2). For small $\mu$ we can approximate the solution of the transcendental equation (34) by assuming $\tanh{\frac{y_{\mu}}{2}}\approx 1$, leading to

[TABLE]

and thus

[TABLE]

Combining Eqs. (27) and (36) we can calculate the limit

[TABLE]

which gives the upper bound on the relative separation between $r^{*}$ and $r^{**}$.

We now show an alternative derivation of the FPT statistics, which should shed some light on the features of the CTRW subdiffusion with stochastic resetting. For a process starting at $t_{0}=0$ from $x_{0}=0$ we can link the distribution of the FAT to a position of the target $x$, denoted in the following as $\rho$, with the free propagator $W$ (without absorbing boundaries) montroll1973random ; shlesinger1974asymptotic ; scher1975anomalous ; chechkin2003first ; kusmierz2015optimal

[TABLE]

In the time-Laplace space (38) simplifies to an algebraic equation

[TABLE]

We can now combine (10) and (39) to obtain the following formula for the Laplace transform of the FAT’s probability density function in the process with hard resets

[TABLE]

with $\rho_{0}(s)$ standing for the (Laplace-transformed) PDF of first arrival times in the (corresponding) process free of resets. Equation (40) has previously been derived Reuveni2016 ; pal2016first using the same general technique that we used here to derive formulas for $\langle T_{r}\rangle$ and $\mathrm{Var}T_{r}$ and one can easily verify that it leads to the same expressions (plug in (8) and compare with (21) and (31)).

IV Second model: resetting of the position

In our second model we assume that resets bring the particle back to $x=0$, but do not affect the waiting times, i.e. after the reset the particle still waits for the next jump as scheduled by the previously generated waiting time. We will refer to this resetting mechanism as soft resets. If the process is space homogeneous (invariant under translations), one can first generate the trajectory of the process without resets, and then introduce the resets by shifting the trajectory from the times of resets onward.

Similarly to the standard CTRW scenario, we derive the formula for the evolution of the PDF of the process $W(x,t)$ by considering a set of renewal integral equations linking the propagator and an auxiliary function $Q(x,t)$, which describes the density of particles jumping from the position $x$ at time $t$. We assume that initially the distribution of the particles is described by $P_{0}(x)\equiv W(x,t=0)$ and all particles have their waiting times generated at $t=0$ (no history). Resets are generated from the exponential distribution, thus the probability that resets do not happen up to time $t$ is given by $\Psi_{R}(t)=\exp(-rt)$.

[TABLE]

where the four terms correspond to the following cases:

1. no jumps and and no resets up to time $t$,
2. no jumps and at least one reset up to time $t$,
3. last jump at time $t-\tau$, no resets from $t-\tau$ up to $t$,
4. last jump at time $t-\tau$, at least one reset from $t-\tau$ up to $t$. Similarly, the equation for $Q(x,t)$ reads

[TABLE]

Notice that the only change on the RHS is the usage of the waiting times PDF $\psi$ instead of the survival probability $\Psi$. We can rewrite (41) and (42) in the Fourier-Laplace space

[TABLE]

[TABLE]

where we have introduced

[TABLE]

which describes the total flow of particles at time $t$ independently from the position. Since resets do not affect waiting times, we can write a simple renewal equation

[TABLE]

which immediately leads to

[TABLE]

Algebraic equations (43), (44), and (47) lead to the following general solution

[TABLE]

As before, we assume that $\psi(s)\approx 1-A_{\mu}s^{\mu}$ and $g(k)\approx 1-D_{0}k^{2}$, take the limits $A_{\mu}\to 0$ and $D_{0}\to 0$ with $A_{\mu}/D_{0}\to D$, and arrive at

[TABLE]

which is our central result in this section.

As we will now show, the process posses rather peculiar qualities. The stationary probability is given simply by

[TABLE]

and is the same as the initial distribution if $x_{0}=0$ (i.e. ${P_{0}(k)=1}$). It does not, however, imply that the process is trivial ($X_{t}=0$). In order to show this we calculate the variance as a function of time. As in the previous example we assume that $x_{0}=0$ and obtain the following expression for the MSD in the Laplace space

[TABLE]

Inverting the Laplace transform leads to the integral expression

[TABLE]

which can be represented in terms of the confluent hypergeometric function of the first kind

[TABLE]

The MSD as a function of time admits maximum value at $t^{*}=z_{0}/r$, where $z_{0}$ is the positive solution of the equation

[TABLE]

which can rewritten in terms of a zero of the confluent hypergeometric function

[TABLE]

The maximum value of the MSD takes the form

[TABLE]

By expanding (51) we see that $\langle X^{2}(s)\rangle\propto s^{-\mu}$ for $s\ll r$, which translates to $\langle X^{2}(t)\rangle\propto 1/t^{1-\mu}$ for large $t$. The stationary state is approached with a slow, algebraic decay, which gets slower with larger values of $\mu$, see Fig. 3.

To sum up, the initial and the stationary distribution are the same, with a nontrivial transient behavior. This bears a resemblance to homoclinic orbits in dynamical systems strogatz2014nonlinear and is related to the stationary state splitting observed recently in a similar model shkilev2017continuous . One may think that this non-monotonic behavior is induced by subdiffusion. In order to show that this is not the case let us modify the model by introducing a power-law, heavy-tailed distribution of jumps, i.e. $g(k)\approx 1-D_{0}|k|^{\alpha}$, with $\alpha\in(0,2)$. The corresponding process without resets is described by the following FFPE

[TABLE]

where $\partial_{|x|}^{\alpha}$ denotes the Riesz fractional space operator metzler2000random . Sample paths of the process are discontinuous, and its increments are described by Lévy $\alpha$-stable distribution fellerBook , hence the name Lévy flights mandelbrotBook ; hughes1981random . Additionally, the variance of the position is infinite for any $\alpha<2$. Nevertheless, one can still differentiate between superdiffusive and subdiffusive cases using quantiles (we could avoid this difficulty by using Lévy walks, but we choose Lévy flights for the sake of simplicity of the calculations). In the case of the process with soft resets, the same calculations as before lead to the following formula

[TABLE]

The underlying process without resets is self-similar and its quantiles scale with time as $\mu/\alpha$. Therefore the process behaves subdiffusively for $2\mu<\alpha$ and superdiffusively for $2\mu>\alpha$. However, it is easy to show that the non-monotonic behavior of the process with position resets appears independently of $\alpha$, since the stationary probability

[TABLE]

which holds for $\mu<1$ and any $\alpha$. This time we cannot use the MSD as a proxy to see the dispersion in time, but we can plot quantile lines to confirm the non-monotonic behavior, see Fig. 4. Thus, it is evident that the non-monotonic behavior of the process relies on the trapping with power-law waiting times. For the same reason no such non-monotonic behavior is observed in the case of subdiffusion with hard resets, wherein power-law waiting times are effectively exponentially truncated.

We hypothesize that this is related to the weak ergodicity breaking, i.e. we expect the non-monotonic behavior will not appear in the fractional Brownian motion or random walks on fractals. If this is the case, one could use the resetting as a computational method to assess whether the process at hand is weakly non-ergodic.

Fig. 4 illustrates another interesting feature of the CTRW with soft resets: Quantile lines converge to zero in a finite time, which means that even for $t<\infty$ there is a finite probability $p_{\mathrm{origin}}(t)$ that the particle is exactly at the origin. This probability corresponds to the singular part of $W(x,t)$ and can be calculated as $\lim_{k\to\infty}W(k,s)$, which for any $g(x)$ that is continuous at $x=0$ reads

[TABLE]

which can be represented in terms of the confluent hypergeometric function

[TABLE]

Fig. 3 shows the behavior of $p_{\mathrm{origin}}(t)$ for different values of $\mu$. Our analytical predictions (61) are fully corroborated by the stochastic simulations. Note that for any value of $\mu<1$ the density is concentrated in the atom in the limit of $t{\to\infty}$. However, the closer the value of $\mu$ to $1$, the slower the convergence, in line with the asymptotic ($t\gg 1/r$) formula

[TABLE]

We now turn our attention to the FAT statistics in the model of subdiffusion with soft resets. Since resets do not affect waiting times for random jumps and for $\mu<1$ these waiting times have infinite expected value, we expect the MFAT to be infinite for $\mu<1$. Moreover, since at any given moment some particles are trapped at the origin and the proportion of trapped particles increases with time, one may expect that there is a finite probability that the target will never be found. This is, surprisingly, not the case. By means of the general formula for the first passage times density (39) together with (49) we find

[TABLE]

which can be rewritten as

[TABLE]

compare to (40). In contrast to hard resets, soft resets do not lead to finite MFAT–as expected $\langle T_{r}\rangle=\infty$ for $\mu<1$. However, $\lim_{t\to\infty}S(t)=0$, i.e. the particle reaches the target almost surely. The tail of $\rho(t)$ behaves like $1/t^{1+\mu}$ and is lighter than the tail of the process without resets $1/t^{1+\mu/2}$.

V External potentials: drift-subdiffusion with resets

We can introduce external potentials into the subdiffusive CTRW and the corresponding FFPE in different ways: trapping (waiting times) may directly affect both random and potential-dependent terms metzler1999anomalous ; metzler2000random . The FFPE takes the form

[TABLE]

where $U(x)$ is the time-independent external potential and ${F(x)\propto-U^{\prime}(x)}$ is the corresponding external force (for a time-dependent generalization see henry2010fractional ). In an alternative approach, waiting times may be affected by the external force, leading to another, more complicated form of the FFPE, in which both terms depend on the external force fedotov2015subdiffusion .

Here we explore (65) with a constant force $f=-U^{\prime}(x)$, i.e. $g(k)\approx 1+ikf_{0}-D_{0}k^{2}$, where $f_{0}\to 0$, $D_{0}\to 0$, $f_{0}/A_{\mu}\to f$, and $D_{0}/A_{\mu}\to D$. The corresponding free propagator (without absorbing boundaries) takes the form

[TABLE]

In order to calculate the FAT statistics we again place an absorbing target at $x=0$ and apply the method of images with a linear combination of solutions (66) of (65) starting from $\pm x_{0}$. For simplicity of the notation and without loss of generality we assume that the process starts from $x_{0}>0$. The formula

[TABLE]

with $c=\sqrt{s^{\mu}/D}$, is a straightforward generalization the known case of diffusion fellerBook . The corresponding survival probability reads

[TABLE]

The infinite-time survival depends on the sign of $f$ and is equal to [math] if $f<0$ (target is found with probability $1$), whereas for positive force the target is found with probability $\exp(-2\bar{f}x_{0})$. The MFAT is finite if and only if $f<0$ and $\mu=1$ and in this case is given by the ballistic time $\langle T_{r}\rangle=x_{0}/|f|$.

The propagator of the corresponding model with hard resets to the origin can be calculated from (10), and the corresponding stationary distribution reads

[TABLE]

where $\bar{f}=\frac{f}{2D}$ and $a=\sqrt{r^{\mu}/D}$. The corresponding MFAT in the process with hard resets to the initial position $x_{0}$ can be calculated with (20) and reads

[TABLE]

It has recently been shown that whether or not resets can lower the MFAT of diffusion with drift is controlled by the Péclet number ray2018p . Our results show that hard resets can always help for $\mu<1$, since in this case the process without resetting has infinite MFAT.

The propagator of the model with soft resets takes the form

[TABLE]

As before, Model II exhibits trivial stationary distribution ${f_{s}(x)=\delta(x)}$ with nontrivial transient. In this case the non-monotonic behavior can be already observed in the first moment

[TABLE]

where the second term behaves like the variance in the model without the drift term (51). One can easily verify that in the case of soft resets drift does not change the main conclusions regarding the FAT statistics: yet again the particle almost surely reaches the target, but the MFAT is infinite.

VI Conclusions

Two models of subdiffusion with stochastic resetting have been presented. Their statistical properties, including stationary distributions, transient behavior of the probability distribution functions, as well as first two moments of first arrival time have been calculated. While the subdiffusion with hard resets provide a rather straightforward generalization of diffusion with stochastic resetting, soft resets lead to quite peculiar features of the process, including non-monotonic behavior of the MSD and non-intuitive properties of the first passage time statistics.

In the following we list a few possible applications of these results.

A number of models of molecular motors dynamics based on ratchet mechanism in a thermally fluctuating environment have been proposed before vale1990protein ; cordova1992dynamics ; astumian1994fluctuation ; reimann2002brownian ; astumian2002brownian ; Klumpp2005 ; RevModPhys.81.387 ; lisowski2015cargo . Since the cytoplasm in the living cells is crowded, the transport inside of them is observed to be subdiffusive (cf. Section I). Therefore, some authors have devised molecular motors models in subdiffusive environments PhysRevLett.85.5655 ; Zaid2009710 ; PhysRevE.80.011912 ; goychuk2014molecular ; goychuk2015anomalous . Our first model of subdiffusion with stochastic resetting may provide an useful building block in such models, wherein resetting would describe binding and unbinding of molecular motors from microtubules or disappearance and appearance of particles due to chemical reactions rotbart2015michaelis ; Reuveni2016 .

The second model introduced in this work may be used as a starting point to construct novel resampling statistical methods efron1982jackknife ; westfall1993resampling ; sauerbrei1999use ; yekutieli1999resampling ; mackinnon2004confidence ; howell2012statistical ; good2013permutation , e.g. for an estimation of the tail index or testing hypothesis whether the observed trajectories are subdiffusive meroz2015toolbox . The idea is to create auxiliary trajectories by randomly resetting positions in the original data, with the hope of the new data being easier to handle, especially when a limited number of trajectories is given. The details of the algorithm and whether this hope is justified will be the subject of further research.

Acknowledgements.

This work was supported by RIKEN Center for Brain Science.