Non-equilibrium steady states of stochastic processes with intermittent resetting

Suppose that you are working on a difficult problem. While trying to find a solution, you may get the impression that you are stuck or got on the wrong track. A natural strategy in such a situation is to reset from time to time and start over. This behavior can be modeled by a stochastic exploration process, which is interrupted by a random resetting to the initial condition. Over the last years a special case of such processes, a diffusion process interrupted at constant rate by reset events, has attracted considerable attention, because it represents a particularly simple and analytically approachable example of a non-equilibrium steady state (NESS) [1–5].

The investigation of processes with random resetting is also one of natural interest to the study of first passage times, e.g. in the catalysis time of chemical reactions [6], in kinetic proofreading [7], and in areas where search optimality is crucial [8]. Furthermore, processes with reset are studied in population dynamics, where resets are interpreted as catastrophic events corresponding to the extinction of the population, followed by a resurgence [9].

All of the previous works, including the extensions to Lévy Flights [10] and fluctuating interfaces with stochastic resetting [11], assume the special case that resets occur at a constant rate γ, i.e. the case in which the waiting times between the resets are exponentially distributed. In this case the resetting procedure is memoryless, which has a straight–forward but interesting implication: the NESS of the process is equivalent to the distribution of end-points of an ensemble of trajectories of the same process without resetting that are trapped at a constant rate γ. Consequently, if the propagator of the process without resetting is known, the NESS can be directly obtained by averaging the propagator over the exponential distribution of inter-reset times [1]. However, the limitation to constant rate resetting severely restricts the applicability to memoryless processes. Even the simple example of a process in which multiple identical steps have to be completed for a reset, which gives rise to Gamma-distributed waiting times, cannot be described within the current framework. Such processes occur e.g. in chemical resetting due to multistep dissociation reactions [12]. Neither can the important case of a rate that depends on the time τ elapsed since the last reset event be captured anymore. Processes of this form can be described using time-dependent rates $\tilde{\gamma }(\tau )=(\alpha \gamma ){(\tau \gamma )}^{\alpha -1}$, which leads to Weibull distributed waiting times. Here the probability for a event decreases (increases) with τ for $\alpha \lt 1$ ($\alpha \gt 1$), thus giving rise to more intermittent (regular) reset sequences, as compared to the constant rate case ($\alpha =1$). The Weibull distribution is widely applied when the conditions for strict history-independence are violated [13].

In this article, we present the theory of stochastic resetting with arbitrary waiting time distributions $\psi (\tau )$, (see figure 1), which includes the aforementioned examples as special cases. While for stochastic processes that reset at a constant rate, the master equation can be immediately formulated, here we first have to derive the governing equation for arbitrary resetting time distributions. The corresponding master equation equation (6) is non-local in time and is one of our main results. Despite the fact that our theory has to include integrals over the complete history of the process, we can explicitly calculate stationary distributions for different waiting time distributions, $\psi (\tau )$, and quantify the temporal relaxation of the moments towards this NESS. We find that the history-dependent processes with a general resetting time distribution have a rich structure and exhibit many new properties that are not present in the case of resetting with a constant rate. Furthermore, we demonstrate that the success of a search depends on the full waiting time distribution and not only on the characteristic resetting time scale as one could expect from the memoryless process. A possible application of our theory is the efficiency of random searchers that are confined by their need to regularly return to a home location. We find that, under this constraint, the search success can be optimized by adapting the distribution of waiting times between the returns.

Zoom In Zoom Out Reset image size

Consider the motion of a particle that, between reset events, is described by the stochastic differential equation

$$\begin{eqnarray}&&\dot{x}(t)=F(x)+\xi (t),\end{eqnarray} \tag{ 1 }$$

where $\xi (t)$ is Gaussian white noise with correlation $\langle \xi (t)\xi ({t}^{\prime })\rangle =2D\delta (t-{t}^{\prime })$. Furthermore, let the waiting time distribution between the resets be $\psi (\tau )$, after which the particle is reset to the location x_r. The probability ${\nu }_{n}(t)$ that the nth reset event happens at time t satisfies the renewal equation ${\nu }_{n}(t)={\displaystyle \int }_{0}^{t}\psi (t-{t}^{\prime }){\nu }_{n-1}({t}^{\prime }){\rm{d}}{t}^{\prime },$ stating that the probability to have the nth event at time t is given by the probability that the $n-1$ th event occurred at time ${t}^{\prime }$ and that the next event occurs after time $t-{t}^{\prime }$ [14]. We assume that the process starts with a reset event. Then, the probability that an event occurs at time t is $\nu (t)=\delta (t)+{\sum }_{n=1}^{\infty }{\nu }_{n}(t)$, where the δ-distribution accounts for the initial condition.

The probability $p(x,t)$ to find the particle at location x at time t is

$$\begin{eqnarray}&&p(x,t)={\displaystyle \int }_{0}^{t}w(t-{t}^{\prime })p(x,{x}_{{\rm{r}}};t-{t}^{\prime })\nu ({t}^{\prime }){\rm{d}}{t}^{\prime },\end{eqnarray} \tag{ 2 }$$

where $p(x,{x}_{{\rm{r}}};t-{t}^{\prime })$ is the transition amplitude of the process defined by equation (1). The distribution $w(t-{t}^{\prime })$ gives the probability that no resetting event occurs between ${t}^{\prime }$ and t and is related to the waiting time distribution $\psi (\tau )$ according to $w(\tau )=1-{\displaystyle \int }_{0}^{\tau }\psi (t){\rm{d}}t$. To derive an evolution equation for $p(x,t)$, we switch to the Laplace domain $\hat{p}(x,s)={\displaystyle \int }_{0}^{\infty }{\rm{d}}t\;{{\rm{e}}}^{-{st}}p(x,t)$. Using the operator representation $p(x,{x}_{{\rm{r}}};t-{t}^{\prime })={{\rm{e}}}^{(t-{t}^{\prime }){ \mathcal L }(x)}\delta (x-{x}_{{\rm{r}}})$ [15] and applying the convolution and shifting theorem, the Laplace transform of equation (2) reads

$$\begin{eqnarray}&&\hat{p}(x,s)=\hat{w}(s-{ \mathcal L }(x))\delta (x-{x}_{{\rm{r}}})\hat{\nu }(s).\end{eqnarray} \tag{ 3 }$$

Noting that the Laplace transform of $\nu (t)$ is

$$\begin{eqnarray}&&\hat{\nu }(s)=1+\hat{\psi }(s)/(1-\hat{\psi }(s)),\end{eqnarray} \tag{ 4 }$$

we can use equations (3) and (4) to obtain

$$\begin{eqnarray}&&s\;\hat{p}(x,s)-{p}_{0}(x)={ \mathcal L }(x)\hat{p}(x,s)+\displaystyle \frac{\hat{\phi }(s)}{s}\delta (x-{x}_{{\rm{r}}})-\hat{\phi }(s-{ \mathcal L }(x))\hat{p}(x,s),\end{eqnarray} \tag{ 5 }$$

where the time evolution kernel $\hat{\phi }(s)$ is given by $\hat{\phi }(s)=s\hat{\psi }(s)/(1-\hat{\psi }(s))$ [16] and ${p}_{0}(x)=\delta (x-{x}_{{\rm{r}}})$ is the initial condition. Finally, Laplace inversion of equation (5) leads to the generalized master equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}p(x,t)={ \mathcal L }(x)p(x,t)+{\displaystyle \int }_{0}^{t}\phi (t-{t}^{\prime })[\delta (x-{x}_{{\rm{r}}})-{{\rm{e}}}^{(t-{t}^{\prime }){ \mathcal L }(x)}p(x,{t}^{\prime })]{\rm{d}}{t}^{\prime }.\end{eqnarray} \tag{ 6 }$$

The first term on the rhs accounts for the dynamics between the resetting events, the second term is the source term and describes the resets to x_r. The third term is a sink term; it describes the density of particles that have propagated to x at time t during the time $t-{t}^{\prime }$ under the influence of the dynamics equation (1) and are then subject to resetting. The integral of the time-evolution kernel corresponds to the time-dependent density of reset events, i.e. $\nu (t)={\displaystyle \int }_{0}^{t}\phi (\tau ){\rm{d}}\tau$. The special case of resetting at a constant rate γ is recovered by choosing an exponential waiting time distribution $\psi (\tau )=\gamma \;{{\rm{e}}}^{-\gamma \tau }$ that leads to $\phi (\tau )=\gamma \delta (\tau )$. Inserting this kernel into equation (6) leads to the Markovian master equation studied in [1]. The solution of equation (6) in Laplace space is

$$\begin{eqnarray}&&\hat{p}(x,s)=\displaystyle \frac{{p}_{0}(x)}{s-{ \mathcal L }(x)+\hat{\phi }(s-{ \mathcal L }(x))}+\displaystyle \frac{\hat{\phi }(s)\delta (x-{x}_{{\rm{r}}})}{s[s-{ \mathcal L }(x)+\hat{\phi }(s-{ \mathcal L }(x))]}.\end{eqnarray} \tag{ 7 }$$

Here the first term on the rhs describes the decay of the initial condition and the second term accounts for the resets. While the Laplace inversion of equation (7) can in general only be carried out by numerical means, we will show that the time-dependent moments of the distribution can be calculated analytically. Remarkably, even the full stationary solution ${p}_{\mathrm{st}}(x)={\mathrm{lim}}_{t\to \infty }p(x,t)$ can be given by employing the Tauberian theorems [17], which allow to consider instead the limit ${\mathrm{lim}}_{s\to 0}\hat{p}(x,s)$ in Laplace space.

3.1. Constant advection with intermittent resetting

Before we come to the diffusion processes with resetting, let us first exemplify equation (6) using the example of constant advection with velocity v between the reset events (figure 1(c)). The discussion of this simplest of all cases helps to illustrate the intricacies that come with the inclusion of non-memoryless resetting. In particular, we set $F(x)=v$ and D = 0 in equation (1) and consider the Liouvillian ${ \mathcal L }(x)=-v{\partial }_{x}$. Without loss of generality let us from now on assume the reset position to be x_r = 0. Using equation (6), we obtain for the time evolution of the first moment in Laplace space

$$\begin{eqnarray}&&\langle \hat{x}(s)\rangle =\displaystyle \frac{v}{s(s+\hat{\phi }(s))}\left(1+\displaystyle \frac{{\rm{d}}}{{\rm{d}}s}\hat{\phi }(s)\right).\end{eqnarray} \tag{ 8 }$$

Clearly, the mean of the process depends on the precise form of the waiting time distribution through $\hat{\phi }(s)$ even in the stationary limit. For concreteness, let us consider the waiting times to be Gamma-distributed with rate γ and shape parameter α, i.e. $\psi (\tau ;\alpha ,\gamma )={\tau }^{\alpha -1}{{\rm{e}}}^{-\gamma \tau }{\gamma }^{\alpha }/{\rm{\Gamma }}(\alpha )$. The mean of this distribution is $\langle \tau \rangle =\alpha /\gamma$ and for $\alpha =1$ it includes the exponential distribution. The shape parameter α regulates how intermittently the resetting events occur. For $\alpha \lt 1$ these events occur more intermittently the smaller α gets, whereas for $\alpha \to \infty$ the resetting events occur regularly with period $\langle \tau \rangle$. In figure 2(a) we compare $\langle x(t)\rangle$ for Gamma-distributed waiting times with the same mean waiting time but with different rate and shape parameters, calculated using equation (8). Interestingly, the asymptotic value of $\langle x(t)\rangle$ depends on the shape of the distribution, even when the mean time between resets, $\langle \tau \rangle$, is kept constant. We can also calculate the asymptotic, steady-state value of the mean by expanding equation (8) in s, keeping only the lowest order, and then performing the Laplace inversion. This leads to ${\langle x\rangle }_{\mathrm{st}}=\displaystyle \frac{v\langle \tau \rangle }{2}\frac{1+\alpha }{\alpha }$, which only for exponentially distributed waiting times is equal to the typical distance ${x}_{s}=v\langle \tau \rangle$ that a particle moves between resets (figure 2(b)). The stationary distribution can be calculated using the same approach. Inserting the Liouvillian ${ \mathcal L }(x)=-v{\partial }_{x}$ in Fourier space, ${ \mathcal L }(k)={\rm{i}}{vk}$, in equation (7), and substituting the appropriate $\hat{\phi }(s)$, the resulting expression can be considered in the limit $s\to 0$. For example, in the case of Gamma-distributed resetting times, the resulting expression in Fourier space is ${p}_{\mathrm{st}}(k)={\rm{i}}\gamma (1-{(1-{\rm{i}}k/\gamma )}^{-\alpha })/(\alpha k)$, which can easily be Fourier inverted. Several examples are given in figure 2(c). The stationary distribution exhibits a transition from a broad-tailed distribution for $\alpha \lt 1$, where rare large excursions dominate the process, to a uniform distribution in the limit $\alpha \to \infty$, where the process is reset deterministically after time $\langle \tau \rangle$.

Zoom In Zoom Out Reset image size

3.2. Diffusion with intermittent resetting

Let us now focus on the case of free diffusion between the resets, i.e. $F(x)=0$, with ${ \mathcal L }(x)=D{\partial }_{x}^{2}$. In this case, the process is symmetric around the reset point and the mean is simply x_r (set to zero for simplicity). Interestingly, calculating the second moment of the process, $\langle {x}^{2}(t)\rangle$ from equation (6), leads to exactly the same equation as the mean of the advection process, but with v replaced by $2D$. Thus, considering the time evolution of the variance and its final state leads to exactly the same curves as shown for the mean of the advection in figures 2(a) and (b), as we have also confirmed numerically. The stationary distribution ${p}_{\mathrm{st}}(x)$ can again be calculated by considering equation (7) in the limit of small s. For the Gamma-distribution, one obtains after Laplace inversion

$$\begin{eqnarray}{p}_{\mathrm{st}}(x) & = & \frac{1}{\sqrt{\pi d\gamma }}{}_{1}{F}_{2}\left(-\frac{1}{2};\frac{1}{2},\frac{1}{2}-\alpha ;\frac{{x}^{2}}{4d\gamma }\right)-\displaystyle \frac{| x| }{2{da}\gamma }\\ & & +\displaystyle \frac{{| x| }^{1+2\alpha }{\rm{\Gamma }}(-2\alpha -1)}{{(d\gamma )}^{\alpha +1}\pi \alpha }{}_{1}{F}_{2}\left(\alpha ;\alpha +1,\alpha +\frac{3}{2};\frac{{x}^{2}}{4d\gamma }\right)\mathrm{sin}(\pi \alpha ),\end{eqnarray} \tag{ 9 }$$

where ${}_{q}{F}_{{\rm{r}}}$ denotes the generalized hypergeometric function. The general solution equation (9) reduces to particularly simple expressions for integer values of α, e.g. $\mathrm{exp}(-| x| /\sqrt{d\gamma })(| x| +3\sqrt{d\gamma })/(8d\gamma )$ for $\alpha =2$. In figures 3(a) and (b) we compare the analytical solution equation(9) to numerical simulations of the process and find a perfect agreement. More intermittent processes, i.e. smaller α lead to broader (but still exponentially decaying) stationary distributions, since longer diffusive excursions become more likely. The stationary standard deviation, ${\sigma }_{\mathrm{st}}=\sqrt{D\langle \tau \rangle (1+\alpha )/\alpha }$ only corresponds to the typical distance diffused by the particle ${x}_{\mathrm{typ}}=\sqrt{2D\langle \tau \rangle }$ for $\alpha =1$. As already mentioned a further important example is the case of drawing the resetting times from a Weibull distribution with density $\psi (\tau ;\alpha ,\gamma )=\gamma \alpha {(\gamma t)}^{\alpha -1}{{\rm{e}}}^{-\gamma {t}^{\alpha }}$, see figures 3(c) and (d). We can again calculate the stationary distribution analytically and observe the same qualitative behavior as for the Gamma distribution. The interesting case of regular resetting at $\langle \tau \rangle$ can be obtained in the limit of $\alpha \to \infty$ of either the Gamma or the Weibull waiting time distribution (while keeping $\langle \tau \rangle$ fixed), or indeed as a limit of any waiting time distribution that approaches $\psi (\tau )=\delta (\tau -\langle \tau \rangle )$ as an appropriate limit. In this important case, where resets are precisely timed, we can express the stationary distribution concisely as

$$\begin{eqnarray}&&{p}_{\mathrm{st}}^{(\mathrm{det})}(\tilde{x})=\frac{1}{\sqrt{D\langle \tau \rangle }}\left[\frac{1}{\sqrt{\pi }}{{\rm{e}}}^{-{\tilde{x}}^{2}}-| \tilde{x}| (1-\mathrm{erf}(| \tilde{x}| )\right],\end{eqnarray} \tag{ 10 }$$

where $\mathrm{erf}$ denotes the error function and $x=2\;\tilde{x}\sqrt{D\langle \tau \rangle }$.

Zoom In Zoom Out Reset image size

Finally, observe that if the resetting is governed by a scale-free waiting time distribution with asymptotic behavior $\psi (\tau )\sim {\tau }^{-1-\delta }$ with $0\lt \delta \lt 1$, the intensity of events decays as $\nu (t)\sim {t}^{\delta -1}$. Therefore, for long times the process is dominated by one large diffusive excursion and does not exhibit a NESS.

Evidently, the resetting protocol also has a strong influence on the search efficiency to find a target located at distance x₀ from the origin. For the case of constant rate resetting, an x₀-dependent optimal rate has been obtained in [1]. However, in many situations the search strategy is constrained by the need to return on average after a characteristic time $\langle \tau \rangle$. For example, many foraging animals have to return home regularly to rest because of exhaustion of energy or to feed their offspring. Our previous results indicate that the efficiency can also be optimized by adapting the functional form of the distribution. In figure 4(a) we show for Gamma-distributed resetting times with a fixed mean resetting time $\langle \tau \rangle$ that the mean first passage time indeed exhibits a x₀-dependent minimum corresponding to an optimal ${\alpha }_{\mathrm{opt}}({x}_{0})$. The existence of such a minimum can be intuitively understood by inspecting the limits of the resetting protocol. For $\alpha \to \infty$, resets occur deterministically after time $\langle \tau \rangle$ and targets at ${x}_{0}\gg {x}_{\mathrm{typ}}$ are almost surely never encountered. On the other hand for $\alpha \to 0$ the reset protocol is extremely intermittent and the process is typically dominated by a single large excursion. In this case the process resembles diffusion without resetting, with a diverging mean first passage time. To determine ${\alpha }_{\mathrm{opt}}({x}_{0})$ we make the reasonable assumption that for ${x}_{0}\geqslant {x}_{\mathrm{typ}}$ the standard deviation of the steady state for ${\alpha }_{\mathrm{opt}}({x}_{0})$ should be proportional to x₀, i.e. ${\sigma }_{\mathrm{st}}[{\alpha }_{\mathrm{opt}}({x}_{0})]/{x}_{0}=C$. We then obtain

$$\begin{eqnarray}&&{\alpha }_{\mathrm{opt}}({x}_{0})={({C}^{2}{x}_{0}^{2}/(D\langle \tau \rangle )-1)}^{-1},\end{eqnarray} \tag{ 11 }$$

where the C has to be determined numerically once and for Gamma-distributed resetting times is $C\approx 0.83$. In figure 4 we compare our analytical prediction equation (11) with the numerically determined values of ${\alpha }_{\mathrm{opt}}({x}_{0})$ and find a very good agreement. For small ${x}_{0}\lt {x}_{\mathrm{typ}}$ equation (11) is no longer valid and it is very hard to identify ${\alpha }_{\mathrm{opt}}$ because the stationary distribution becomes increasingly insensitive to variation of α for $\alpha \gt 1$, see figure 4(a). These results indicate that an intermittent return strategy should be adopted when targets remote from the home may exist. While, we have focused here on optimal search in one dimension, extending our work to two and higher dimensions poses no new conceptual difficulties, and it would be interesting to compare such an extension to recent results on intermittent search strategies [18].

Zoom In Zoom Out Reset image size

In conclusion, our study reveals that the NESS as well as the search efficiency depend on the precise form of the resetting distribution. Interestingly, for continuous time random walks—a related class of models, where the times between successive jumps are drawn from a waiting time distribution—such a dependence of the long time behavior on the precise form of the waiting time distribution is not observed. Here the knowledge of the characteristic waiting time is sufficient to characterize the steady state distribution of the process in a confining potential [19]. Finally, we note that, for general resetting times, the steady state of the process cannot be simply obtained by averaging the propagator over the waiting time distribution. Instead, one has to average over the distribution of times that have elapsed since the last reset event, the so-called backward occurrence time distribution. This distribution is only identical to the waiting time distribution for memoryless processes, i.e. exponentially distributed waiting times. For other waiting time distributions the corresponding backward occurrence time distribution can in general not be obtained in a closed form. Our work opens several further interesting directions for future research. Our theory can be used to consider deterministic dynamics with resetting at random locations, and thus has direct relevance for spiking dynamics in neurons that emit spikes upon reaching a random threshold. Furthermore, it would be interesting to investigate the search efficiency of diffusion processes with generalized resetting in the presence of a binding potential.