Exact enumeration approach to first-passage time distribution of non-Markov random walks

TL;DR

This paper introduces an exact enumeration method to analytically derive the first-passage time distribution for non-Markov random walks, including fractional Brownian motion, applicable to various scientific fields.

Contribution

It presents a novel exact enumeration approach for calculating FPT distributions in non-Markov processes, extending analytical tools beyond Markovian assumptions.

Findings

Derived exact FPT distribution for non-Markov walks

Applied method to fractional Brownian motion with Hurst exponent H

Demonstrated applicability to physics, biology, and geology phenomena

Abstract

We propose an analytical approach to study non-Markov random walks by employing an exact enumeration method. Using the method, we derive an exact expansion for the first-passage time (FPT) distribution for any continuous, differentiable non-Markov random walk with Gaussian or non-Gaussian multivariate distribution. As an example, we study the FPT distribution of a fractional Brownian motion with a Hurst exponent $H\in(1/2,1)$ that describes numerous non-Markov stochastic phenomena in physics, biology and geology, and for which the limit $H=1/2$ represents a Markov process.

Figures

Equations

$$\frac{\partial x(t)}{\partial t}=v(t),C(t,t^{\prime })=⟨v(t)v(t^{\prime })⟩.$$

$${\int }_{x_c-vdt}^{x_c}P(x,v,t|x_0,v_0,t_0)\mathrm{d}x=vP(x_c,v,t|x_0,v_0,t_0)\mathrm{d}t,$$

$$n_1(x_c,t_1|x_0,v_0,t_0)={\int }_0^{\infty }vP(x_c,v,t_1|x_0,v_0,t_0)\mathrm{d}v.$$

$$n_p(x_c,t_p;\cdots ;x_c,t_1|x_0,v_0,t_0)={\int }_0^{\infty }\mathrm{d}{v}_p\cdots {\int }_0^{\infty }\mathrm{d}{v}_1{v}_p\cdots v_{1}P(x_c,v_p,t_p;\cdots ;x_c,v_1,t_1|x_0,v_0,t_0).$$

$$\Big{|}\bigcap_{i=1}^{n}\bar{A}_{i}\Big{|}=\Big{|}U-\bigcup_{i=1}^{n}A_{i}\Big{|}=|U|-\sum_{i=1}^{n}|A_{i}|+\sum_{1\leq i\leq j\leq n}|A_{i}\cap A_{j}|-\dots+\sum_{1\leq i\leq j\leq\cdots\leq n}(-1)^{n-1}|A_{i}\cap A_{j}\cap\dots\cap A_{n}|\;,$$

$$f(t)$$

$$=$$

$$f(t)={\psi }^{\prime }(t)e^{-\psi (t)}.$$

$$f(t)\approx n_1(t)\exp [-{\int }_0^t{n}_1(t^{\prime })\mathrm{d}{t}^{\prime }].$$

$${\psi }_{\mathrm{Str}}(t)=-{\int }_0^t{n}_1(t^{\prime })\frac{\ln [1-{\int }_0^{t}R(t,t^{\prime })n_1(t^{\prime })\mathrm{d}{t}^{\prime }]}{{\int }_0^{t}R(t,t^{\prime })n_1(t^{\prime })\mathrm{d}{t}^{\prime }}\mathrm{d}{t}^{\prime }$$

$$⟨x_H(t_1)x_H(t_2)⟩=\frac{1}{2}(|t_1{|}^{2H}+|t_2{|}^{2H}-|t_2-t_1{|}^{2H}),$$

$$f(t)\sim t^{H-2}.$$

$$n_1(t)=p(x_c){\int }_0^{\infty }vp(v|x_c)\mathrm{d}v,$$

$$n_1(t)=\frac{{\Gamma }^2}{2\pi H{t}^{2H-1}}\exp \left(-\frac{y^2}{2}\right)\left\{\frac{H^2{t}^{2H-2}}{{\Gamma }^2}\exp \left(-\frac{y^2{\Gamma }^2}{2}\right)+\frac{H^2{x}_c}{2\Gamma }{t}^{H-2}\sqrt{2\pi }\left[1+\mathrm{erf}\left(\frac{y\Gamma }{\sqrt{2}}\right)\right]\right\},$$

$$n_2(t_1,t_2)$$

$$=$$

$$⟨x_c^{\prime }|x_c,v⟩=x_c\frac{⟨x^{\prime }x⟩}{t^{2H}}+(v-⟨v|x_c⟩)\frac{⟨x^{\prime }v⟩-⟨x^{\prime }x⟩/2t^{2H}}{s_{v|x_c}}.$$

$$C_H(\tau ,\delta )=\frac{{\sigma }^2{\delta }^{2H-2}}{2}\left[{\left(\frac{|\tau |}{\delta }+1\right)}^{2H}+{\left|\frac{|\tau |}{\delta }-1\right|}^{2H}-2{\left|\frac{\tau }{\delta }\right|}^{2H}\right],$$

$$B_H(t)=B_H(0)+\frac{1}{\Gamma (H+1/2)}\left\{{\int }_{-\infty }^0\left[(t-s{)}^{H-1/2}-(-s{)}^{H-1/2}\right]\mathrm{d}W(s)+{\int }_0^t(t-s{)}^{H-1/2}\mathrm{d}W(s)\right\},$$

$$G_H(t)$$

$$+$$

$${\mathbf{z}}^{\prime }\mathbf{Cz}=\sum _{i=1}^n\sum _{i=j}^n{z}_i{z}_j{C}_{i,j}\ge 0\forall z_1,\cdots z_n\in \mathrm{I}\mathrm{R},$$

$$C(t_1,t_2)=\frac{1}{2}(|t_1{|}^{2H}+|t_2{|}^{2H}-|t_2-t_1{|}^{2H}),$$

$$\Phi (t_2,r)=(t_2-r{)}^{{\alpha }_{+}-1/2}-(-r{)}^{{\alpha }_{+}-1/2},$$

$${\int }_{-\infty }^{\infty }\Phi (t_2,r)\Phi (t_1,r)\mathrm{d}r=\kappa C(t_1,t_2)\forall (t_1,t_2)\ge 0,$$

$$\sum _{i=1}^n\sum _{j=1}^n{z}_i{z}_j{C}_{t_i,t_j}$$

$$P_n(\mathbf{Q})=\frac{1}{(2\pi {)}^{n/2}\sqrt{\det {\hat{\mathbf{C}}}_n}}\exp \left(-\frac{\mathbf{Q}{\hat{\mathbf{C}}}_n^{-1}\mathbf{Q}}{2}\right).$$

$$n_1(t)=\frac{{\Gamma }^2}{2\pi H{t}^{2H-1}}\exp \left(-\frac{y^2}{2}\right)\left\{\frac{H^2{t}^{2H-2}}{{\Gamma }^2}\exp \left(-\frac{y^2{\Gamma }^2}{2}\right)+\frac{H^2{x}_c}{2\Gamma }{t}^{H-2}\sqrt{2\pi }\left[1+\mathrm{erf}\left(\frac{y\Gamma }{\sqrt{2}}\right)\right]\right\},$$

$$n_2(t,t^{\prime })={\int }_0^{\infty }v\mathrm{d}v{\int }_0^{\infty }\mathrm{d}{v}^{\prime }{v}^{\prime }p(x_c,x_c^{\prime },v,v^{\prime })=p(x_c){\int }_0^{\infty }\mathrm{d}vvp(v|x_c)p(x_c^{\prime }|x_c,v){\int }_0^{\infty }\mathrm{d}{v}^{\prime }{v}^{\prime }p(v^{\prime }|x_c^{\prime },x_c,v).$$

$$⟨x^{\prime }x⟩$$

$$⟨v^{\prime }x⟩$$

$$⟨x^{\prime }v⟩$$

$$⟨v^{\prime }v⟩$$

Full text

Exact enumeration approach to first-passage time distribution of non-Markov random walks

Shant Baghram,1, Farnik Nikakhtar,1,2, M. Reza Rahimi Tabar,1,3,† Sohrab Rahvar,1,4 Ravi K. Sheth,2 Klaus Lehnertz,5,6,7 and Muhammad Sahimi8,‡

1*Department of Physics, Sharif University of Technology, Tehran 11155-9161, Iran

2Center for Particle Cosmology, University of Pennsylvania, 209 S. 33rd Street, Philadelphia, PA 19104, USA

3Institute of Physics, Carl von Ossietzky University of Oldenburg, Carl von Ossietzky Straße 9-11, 26111 Oldenburg, Germany

4 Department of Physics, College of Science, Sultan Qaboos University, P.O. Box 36, P.C. 123, Muscat, Sultanate of Oman

5Department of Epileptology, University of Bonn, Sigmund Freud Straße 25, 53105 Bonn, Germany

6Helmholtz Institute for Radiation and Nuclear Physics, University of Bonn, Nussallee 14-16, 53115 Bonn, Germany

7Interdisciplinary Center for Complex Systems, University of Bonn, Brühler Straße 7, 53175 Bonn, Germany

8Mork Family Department of Chemical Engineering and Materials Science, University of Southern California, Los Angeles, California 90089-1211, USA*

We propose an analytical approach to study non-Markov random walks by employing an exact enumeration method. Using the method, we derive an exact expansion for the first-passage time (FPT) distribution for any continuous, differentiable non-Markov random walk with Gaussian or non-Gaussian multivariate distribution. As an example, we study the FPT distribution of a fractional Brownian motion with a Hurst exponent $H\in(1/2,1)$ that describes numerous non-Markov stochastic phenomena in physics, biology and geology, and for which the limit $H=1/2$ represents a Markov process.

I. INTRODUCTION

The concept of first passage refers to the crossing of a prespecified location, or some sort of a threshold, in a stochastic trajectory [1]. The distribution of the first-passage times (FPTs), which represents the probability of crossing the trajectory at a specific time or location [2,3] and depends on the nature of the stochastic process, plays a fundamental role in the theory of stochastic processes, as well as in their applications. The FPT distribution makes it possible to investigate quantitatively the uncertainty in the properties of a stochastic system within a finite time. Two important applications are the extinction time of a disease in the models of epidemic phenomena, and the time for a species to reach a critical threshold in population dynamics. In addition, the statistics of the FPT distribution have many applications to diffusion-limited processes in physics [1], chemistry [4], and biology [5], spreading of electrical blackouts [6], epidemiology [7], and even foraging animals [8,9], as well as to understanding transport processes in disordered materials [10], porous media [11,12], neuroscience [13-16], spreading of computer viruses [17], target search processes [18], economics [19], mathematical finance [20,21], psychology [22], cosmology [23,24], and the reliability theory [25]. Through a suitable boundary the FPT presents the first time that the error in the so-called clock model [26] becomes too large and uncontrollable. Rapid detection of anomalies is closely related to recognizing the optimal stopping time of a diffusion process [27] and, hence, the FPT distribution.

Due to their very large number of applications, the FPT properties have been studied extensively, and are well understood when the stochastic phenomena represent a Markov process. As a general rule, however, the dynamics of a given stochastic process in complex media is the result of its interactions with the environment around it, which may contain trapping sites, obstacles, moving parts, active pumps, etc. [28], and cannot be described as a Markov process. Indeed, although the evolution of the set of all microscopic degrees of freedom of a system is Markovian, the dynamics restricted only to the random walker is not [3,29,30]. Experimental realizations of non-Markov dynamics include diffusion of tracers in crowded narrow channels [31] and in complex fluids, such as nematics [32] and viscoelastic solutions [33,34], as well as the dark matter halo mass function [35]. Even in simple fluids, hydrodynamic memory influences various phenomena and, thus, non-Markov dynamics has been reported recently [36].

Using inclusion-exclusion principle and an exact enumeration method, we derive in this paper the FPT distribution of a non-Markov random walk by assuming that the trajectory of the walk is differentiable at every point. As an example, we derive the FPT distributions of fractional Brownian motion (FBM) with a given Hurst exponent $H\in(0.5,1)$. The analytical results are confirmed by extensive numerical simulation and the analysis of $10^{6}$ trajectories.

The rest of this paper is organized as follows. In the next section we describe the exact enumeration approach to derive the FPT probability density of a non-Markovian random walk. We then drive in Sec. III an analytical expression for the FPT distribution of FBM. The results of numerical simulations are presented in Sec. IV, while the paper is summarized in Sec. V. In the Appendix, we provide the details of the derivation of our results.

II. EXACT ENUMERATION METHOD FOR THE FPT PROBABILITY DENSITY

We define a general dynamical equation for a random walk, $x(t)$, driven by a correlated, nonstationary noise (velocity) $v(t)$,

[TABLE]

$x(t)$ is assumed to be continuous and its derivative (velocity) $v(t)$ to be well-defined at any time [37]. The noise $v(t)$ has a zero mean and an arbitrary $n$-point joint distribution $p(v_{n},t_{n};\cdots v_{1},t_{1};v_{0},t_{0})$. The correlation function $C(t,t^{\prime})$ depends on both $t$ and $t^{\prime}$. Because $x(t)$ is a stochastic process, each of its realizations reaches a given barrier $x=x_{c}$ for the first time at a different time $t$, giving rise to a FPT probability density $f(t)$. Consider the trajectories with the initial conditions $x(t_{0})=x_{0}$ and $\dot{x}(t_{0})=v(t_{0})=v_{0}$, crossing the barrier $x_{c}$ in the time interval $t$ and $t+{\mathrm{d}}t$ with $v(t)>0$. The crossing is equivalent to the conditions that $x(t)<x_{c}$ and $x_{c}<x(t+{\mathrm{d}}t)$ [16,38]. If $x_{c}$ is constant, $x(t)$ will lie in the interval $x_{c}-v{\mathrm{d}}t<x(t)<x_{c}$. Then, the probability that $x(t)$ satisfies the passage condition $x_{c}-v{\mathrm{d}}t<x(t)<x_{c}$ is

[TABLE]

where we kept the terms up to the order of ${\mathrm{d}}t$. Since $v(t)>0$ at $x_{c}$, we should integrate over all positive velocities. Therefore, the probability of crossing the barrier $x_{c}$ per unit time is given by [16],

[TABLE]

Equation (2) represents the rate of up-crossing, rather than a density function and, thus, it is not normalized. We generalize Eq. (2) to the joint probability of multiple up-crossings, i.e., $x(t)$ crossing the barrier in each of the intervals $(t_{1},t_{1}+{\mathrm{d}}t),\cdots,(t_{p},t_{p}+{\mathrm{d}}t)$, by integrating over all the crossing points $t_{1},t_{2},\cdots,t_{p}$,

[TABLE]

Using Bayes’ theorem, one may substitute the conditional probability density in Eq. (3) with the joint probability density. In Fig. 1 typical trajectories, as well as the FPT distribution of the FBM for $x_{c}=1$ with $x_{0}=0$ are presented. The trajectories are constructed using the Cholesky decomposition (see below).

A trajectory can cross $x_{c}$ several times (see the lower panel of Fig. 1). We relate the FPT distribution to the statistical properties of the up-crossings, which are considered as point processes with rates $n_{p}$, where $p$ refers to the number of up-crossing. To this end, we look for the fraction of all the trajectories that up-cross $x_{c}$ for the first time at time $t$ with the initial conditions $(x_{0},v_{0})$ at time $t_{0}$, and enumerate them in terms of $n_{p}$. To simplify the notation, we drop $x_{c}$ and the initial conditions.

The rate $n_{1}(t)$ is over-counted through the trajectories that had an up-crossing at shorter times $t_{1}<t$. Therefore, we subtract their fraction from the first term. This stems from the fact that $n_{1}(t)$ is a local function in $t$, but there is no guaranty that a trajectory has not up-crossed before $t$. The over-counting implies that the main problem is a combinatorial counting. Thus, as an enumeration technique we use the inclusion-exclusion principle, one of the most useful principles of counting in combinatorics and probability. According to De Morgan’s laws, in the general and complementary form, the principle of inclusion-exclusion for finite sets $A_{1},\;A_{2},\dots,\;A_{n}$ is expressed by

[TABLE]

where $U$ is a finite universal set containing all the $A_{i}$, and $\bar{A_{i}}$ are the complement of $A_{i}$ in $U$. That the trajectories cross $x_{c}$ for the first time at time $t$ implies that they should not have been crossed at $x_{c}$ at shorter times. We consider $n_{1}(t)$ as the universal set, and define the next subset by $A_{i}=n_{2}(t,t_{i})$, denoting the fraction of trajectories for which the up-crossing at time $t$ is not for the first time, and that they had a previous up-crossing at a shorter time $t_{i}<t$. Then, the FPT distribution is given by, $|\bigcap_{i=1}^{n}\bar{A}_{i}|$, because only the trajectories that have a first up-crossing at time $t$ and do not belong to the subsets $A_{i}$ are of interest. Using Eq. (4), we obtain the FPT distribution [16]:

[TABLE]

where $n_{p+1}(t,t_{p},\dots,t_{1})$ are given by the conditional probabilities (3). The factor $1/p!$ accounts for the number of permutations of the variables $t_{p},\dots,t_{1}$, with the signs explained in Table I. To calculate $n_{p}(t_{p},\dots,t_{1})$, we consider the trajectories in the absence of $x_{c}$ and let them return after an up-crossing and, then, up-cross the barrier $p$ times. The correct counting of such multiple crossings yields the distribution $f(t)$ of the FPT. Equation (5) provides us with the exact expansion of the FPT distribution for any continuous, differentiable non-Markov random walk with Gaussian or non-Gaussian multivariate distribution [16]. We note that a naive truncation of the series would give rise to a non-normalized (diverging in the long-time limit) distribution [39].

Let us define as a point process the time scales at which the trajectories cross $x_{c}$. The distributions of such a point process are the aforementioned rate functions. Since the trajectories have nonzero velocities, successive up-crossings cannot be too close, so that $n_{p}(t_{p},\dots,t_{1})$ is zero if two of its arguments are equal. Such point processes represent systems of nonapproaching random points [40]. There are two types of decoupling approximations to deal with the infinite series in Eq. (5), which are based on approximating the higher-order terms by the lower-order ones and are known as the Hertz and Stratonovich approximations. The general expression for $f(t)$ is given by

[TABLE]

The Hertz approximation is based on assuming that all the up-crossings are independent of each other, and that the correlations between them are negligible. This leads to the following FPT distribution with $\psi_{\rm Hertz}(t)=\int_{0}^{t}n_{1}(t^{\prime}){\mathrm{d}}t^{\prime}$ [24,39,41]:

[TABLE]

In the Hertz approximation $n_{p}(t_{p},\dots,t_{1})$ factorizes to $n_{1}(t_{p})\dots,n_{1}(t_{1})$. In the Stratonovich approximation, we calculate exactly the first and the second terms of the expansion and approximate all the higher-order terms by the first two [35], with the corresponding FPT distribution being in the form of Eq. (6) with [39],

[TABLE]

where $R(t_{i},t_{j})=1-n_{2}(t_{i},t_{j})/[n_{1}(t_{i})n_{1}(t_{j})]$. For simplicity and in order to derive an expression for $f(t)$, we assume in the following that the velocity distribution is Gaussian.

III. ANALYTICAL DERIVATION OF THE FPT DISTRIBUTION OF THE FBM

We now derive the FPT distribution of the FBM with a Hurst exponent $H\in(0.5,1)$, which is defined in terms of its nonstationary correlation function [42]:

[TABLE]

which is positive semidefinite (see the Appendix) with its first derivative (velocity) being the fractional Gaussian noise (FGN) $v_{H}(t)$ so that, $\dot{x}_{H}(t)=v_{H}(t)$. Using physical arguments [43,44], as well as rigorous analysis [45], it was shown that the scaling behavior of the FPT distribution of a FBM has the following long-time behavior

[TABLE]

Given that the FBM and FGN have Gaussian distributions for $x$ and $v$, respectively, we determine $n_{1}(t)$ and $n_{2}(t_{1},t_{2})$ and, therefore, $R(t_{1},t_{2})$ and the FPT distribution in the Hertz and Stratonovich approximations. It is straightforward to show that $n_{1}(t)$ is given by the following expression,

[TABLE]

where $p(v|x_{c})$ is a Gaussian distribution with mean $\langle v|x_{c}\rangle=x_{c}\langle vx\rangle/\langle x^{2}\rangle=x_{c}H/t$ and variance $s_{v|x_{c}}=H^{2}t^{2H-2}/\Gamma^{2}$, where $\Gamma^{2}=\gamma^{2}/(1-\gamma^{2})$ and $\gamma^{2}=\langle xv\rangle^{2}/\langle x^{2}\rangle\langle v^{2}\rangle$. For the FBM, $\langle xv\rangle^{2}=H^{2}t^{4H-2}$ and $\langle x^{2}\rangle=t^{2H}$. In the Appendix, we present an expression for $\langle v^{2}\rangle$ in terms of the Hurst exponent $H$. We find that the explicit expression for $n_{1}(t)$ is given by

[TABLE]

where $y=x_{c}/t^{H}$. Using Eq. (7) we obtain the FPT distribution in the Hertz approximation, which, in general, is accurate for estimating the first peak of the FPT distribution, but it over- or underestimates its tail. Similarly, we find that,

[TABLE]

where all the distributions in Eq. (12) are Gaussian. For example, $p(x_{c}^{\prime}|x_{c},v)$ has the mean (see the Appendix for the variance)

[TABLE]

The correlation functions $\langle x^{\prime}x\rangle$ and $\langle x^{\prime}v\rangle$ are given by Eq. (A.9) in the Appendix, and $s_{v|x_{c}}=H^{2}t^{2H-2}/\Gamma^{2}$. Having $n_{1}(t)$ and $n_{2}(t_{1},t_{2})$ enables one to determine $R(t_{1},t_{2})$ and $f(t)$ in the Hertz and Stratonovich approximations.

IV. NUMERICAL RESULTS

Figure 1 presents the trajectories of a FBM process using the Cholesky decomposition [46,47] (see the Appendix) and their FPT distribution for $H=0.8$, $x_{c}=1$, and $x_{0}=0$. In Fig. 2 the FPT distributions of the FBM trajectories is plotted. The FPT is obtained from the Cholesky method. In these plots, we also show the FPT distributions in the Hertz approximation, which deviate from the FPT directly computed using trajectories. For comparison, the theoretically-predicted tails of the distributions, i.e., $f(t)\sim t^{H-2}$, are also plotted. Figure 2 indicates that the theoretical tails of the FPT in the long-time limit coincide with the FPT distributions computed using the trajectories. As already mentioned above the Hertz approximation predicts correctly the location of the peak of the FPT distribution, but underestimates the tails.

To derive the FTP distribution in the Stratonovich approximation with $H=0.8$, one must calculate $\psi_{\rm Str}(t)$ via Eq. (9), and then use Eq. (6). To avoid any error from the numerical differentiation of $\psi_{\rm Str}(t)$, we determine the integrated FPT distribution via the term $\exp[-\psi_{\rm Str}(t)]$. In Fig. 3 the cumulative FPT distribution is presented for $H=0.6$ and $H=0.8$, indicating that the Hertz approximation deviates clearly from the results computed via the Cholesky decomposition. As shown in Fig. 2, the tail of $f(t)$ in the Hertz approximation does not coincide completely with those obtained by the Cholesky decomposition. Higher-order approximations, e.g., the Stratonovich approximation, are therefore needed, implying that $n_{2}(t_{i},t_{j})$ should not be factorized as $n_{1}(t_{1})n_{1}(t_{2})$. As shown in Fig. 3, the Stratonovich approximation provides better estimations for the FPT distributions.

One may define various measures to study the interdependence of the up-crossing events. The simplest measure is the Fano factor. Consider a time window $T$ and count the mean number (and its variance) of up-crossing events for trajectories in the window. The Fano factor ${\cal F}(T)$ is defined as the variance of the number of up-crossing events in $T$, divided by its mean number, and is written in terms of $n_{2}(t_{1},t_{2})$ and $n_{1}(t)$ [48]. More specifically, the Fano factor is given by ${\cal F}=\langle\Delta N^{2}\rangle/\langle N\rangle$ (with $\langle\Delta N^{2}\rangle=\langle N^{2}\rangle-\langle N\rangle^{2}$), where $\langle N\rangle=\int_{0}^{T}n_{1}(t){\mathrm{d}}t$ and $\langle N^{2}\rangle=\langle N\rangle+\int_{0}^{T}\int_{0}^{T}n_{2}(t_{2},t_{1}){\mathrm{d}}t_{2}{\mathrm{d}}t_{1}$ [48]. For independent point processes, i.e., $n_{2}(t_{2},t_{1})=n_{1}(t_{2})n_{1}(t_{1})$, one has ${\cal F}=1$. Therefore, for a Poisson process ${\cal F}(T)=1$. By definition, ${\cal F}(T)>1$ and ${\cal F}(T)<1$ refer, respectively, to over- and under-dispersion [49]. We plot in Fig. 4 the Fano factor versus the size of the time window $T$, which indicates that, in the long-time limit, the up-crossing point processes are strongly over-dispersed. This means that in such time scales $n_{2}(t_{i},t_{j})$ should not be factorized as $n_{1}(t_{1})n_{1}(t_{2})$, and that the Hertz approximation is not appropriate for estimating the tails of the FPT distribution. A very crucial point to indicate is that in the time span which ${\cal F}\sim 1$ the Hertz approximation works well and it is very near to the Cholesky-derived FPT distribution. On the other and if the Fano factor deviates from unity, it is certain that the Hertz approximation is not suitable for FPT, however this parameter can not quantify the accuracy of Stratonovich approximation.

V. SUMMARY

Except for the limiting case of Markov processes, no exact analytical expression for the FPT distribution of general non-Markov random walks had been derived. In principle, the FPT distribution of non-Markov processes may be obtained from the solution of the associated Fokker-Planck equation with absorbing boundaries in higher dimensions, resulting from the Markovian embedding of a non-Markov process [50]. Even the calculation of the mean FPT for a non-Markov process is, however, a rather difficult task, since the corresponding boundary problem cannot be treated in a straightforward manner [51-55]. We presented a general method for deriving such analytical expressions for the FPT distribution. This is done by using an exact enumeration method based on combinatorics and the inclusion-exclusion principle, which can be generalized to include the FTP distribution of non-Markov random walks in higher dimensions. As an example, analytical results were presented for the FBM with the Hurst exponent $H\in(0.5,1)$, which is a non-Markov process with infinite-range memory, and has wide applications in many disciplines [28]. The numerical results were also compared with two well-known approximations, namely, the Hertz and Stratonovich approximations, which revealed their shortcomings.

ACKNOWLEDGMENTS

S.B. acknowledges the partial support of Sharif University of Technology, Grant No. G960202 for this paper. F.N. acknowledges support from the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1845298.

APPENDIX

We provide the details of the main results presented in the main text of the paper.

A. Variance of the velocity of fractional Brownian motion

The time derivative (increments) of the FBM is the FGN, and has the following correlation function

[TABLE]

where $0<H<1$, and $\tau=t_{2}-t_{1}$. Here, $\delta>0$ is used for smoothing the FBM to make it numerically differentiable [42]. We note that in the limit $\tau\to 0$, the $\delta$-dependence of $\gamma^{2}=\langle xv\rangle^{2}/\langle x^{2}\rangle\langle v^{2}\rangle$ drops out. In the literature [42], there is no unique expression for $\langle v^{2}\rangle$. Here, by generating the FBM trajectories and numerically differentiating them for $H\in(0.5,1)$, the best fit is found to be $\langle v^{2}\rangle=c_{0}+c_{1}H^{m}$, where $c_{0}=-2.47\pm 0.01$, $c_{1}=2.88\pm 0.05$, and $m=-4.72\pm 0.02$.

B. Fractional Gaussian noise

The stochastic representation of the FBM is given by,

[TABLE]

where $dW(s)$ is a Wiener process that is written in terms of the Gaussian white noise $\xi(s)$ as, $dW(s)=\xi(s)ds$. The FGN is then defined by, $G_{H}(t)=dB_{H}(t)/dt$. Taking the time derivative of Eq. (A2) yields

[TABLE]

The second term on the r.h.s. of Eq. (A3) is not finite for $H\in(0,0.5)$. Therefore, the FBM has no well-defined ”velocity” for the Hurst exponent in the range of $(0,0.5)$.

C. Proof for the variance of the FBM being positive semidefinite

A symmetric $n\times n$ real matrix C is the covariance of some random (Gaussian) vector, if and only if it is positive semidefinite, which means that

[TABLE]

where $z$ here is the aforementioned random vector. The FBM has a vanishing mean ($x(0)=0$), while its covariance is given by Eq. (10) of the main text, for $(t_{1},t_{2})\geq 0$ and $H\in(0,1)$. We show that

[TABLE]

is a covariance function. Consider the function

[TABLE]

defined for all $t_{2}\geq 0$ and $r\in{\rm I\!R}$, where $\alpha_{+}={\rm max}(0,H)$ for all $H\in{\rm I\!R}$. Since $H<1$, we can determine $\int_{-\infty}^{\infty}|\Phi(t_{2},r)|^{2}{\mathrm{d}}r<\infty$ and

[TABLE]

where $\kappa$ is a positive and finite constant that depends only on $H$. Therefore, we find

[TABLE]

D. Analytical expressions for ${\bf n_{1}(t)}$ and ${\bf n_{2}(t,t^{\prime})}$ with a Gaussian velocity

Due to the linearity of the system, all the joint probability densities are Gaussian and have the form

[TABLE]

Here, ${\bf Q}=(q_{1}(t_{1}),\dots,q_{n}(t_{n}))$ is an $n$-dimensional vector whose $i$th component is the coordinate $x(t_{i})$ or the velocity $v(t_{i})$ at time $t_{i}$, and $\hat{{\bf C}}_{n}$ is the symmetric $n\times n$ correlation matrix whose entries are the correlation functions between the corresponding components of the vector ${\bf Q}$: $C_{ij}=C_{ji}=\langle q_{i}(t_{i})q_{j}(t_{j})\rangle$. Then, $n_{1}(t)$ is obtained in closed analytical form:

[TABLE]

where $y=x_{c}/t^{H}$ and $\Gamma^{2}=\gamma^{2}/(1-\gamma^{2})$. For the joint densities of multiple up-crossings $n_{p}(t_{p},\dots,t_{1})$ no closed expression can be obtained. We evaluate the integral over $v_{1}$ in Eq. (3) analytically and then perform numerical integration of the resulting expression over $v_{2},\dots,v_{p}$ to determine $n_{p}(t_{p},\dots,t_{1})$. The integrals over time in the expressions for $f(t)$ are also evaluated numerically. For $n_{2}(t,t^{\prime})$, we compute the mean and variance of the conditional distributions,

[TABLE]

Assuming that $t^{\prime}>t$, the correlations are given by

[TABLE]

where, for example, $\langle x^{\prime}x\rangle=\langle x(t^{\prime})x(t)\rangle$. We also know that

[TABLE]

Note that all the conditional distributions are Gaussians and, therefore, they are specified by their mean and variance. For example, for $p(x^{\prime}_{c}|x_{c},v)$ we have

[TABLE]

[TABLE]

For $p(v^{\prime}|x^{\prime},x,v)$, we should calculate the mean and variance of $p(v^{\prime}|x,v)$, which are given by

[TABLE]

[TABLE]

Now, for $p(v^{\prime}|x^{\prime},x,v)$ we obtain

[TABLE]

[TABLE]

E. The Cholesky decomposition

To compute the non-Markovian first up-crossing distribution for the FBM, we must generate trajectories with the correct ensemble properties. Here, we describe an algorithm to generate such trajectories. Equation (9) of the main text defined, $C_{ij}\equiv C(t_{i},t_{j})=\langle x(t_{i})x(t_{j})\rangle$, the correlation between the $x(t)$ between times $t_{i}$ and $t_{j}$. The matrix C is real, symmetric, and positive-definite and, therefore, it has a unique decomposition, ${\bf C}={\bf LL}^{\rm T}$, in which L is a lower triangular matrix, which is known as the Cholesky’s decomposition. We use L to generate the ensemble of the trajectories as follows.

First, consider a vector $\mathbf{\xi}$, which is Gaussian white noise with zero mean and unit variance (i.e. $\langle\xi_{m}\xi_{n}\rangle=\delta_{mn}$). If we generate the desired trajectories as

[TABLE]

then $x_{i}$ will have the correlations of random walk given by

[TABLE]

Since L is triangular, only a sum over $j\leq i$ is needed in matrix calculations, so the method is fast. Next, we provide a proof of the Cholesky decomposition, and present it in terms of the correlation matrix. In a more general context, there is a sufficient condition for a square matrix to have a LU decomposition, ${\bf C}={\bf LU}$, where L and U are, respectively, the lower and upper triangular matrices of C. If all the $n$ leading principal minors of the $n\times n$ matrix C are nonsingular, then C has an LU decomposition. Let us recall that the $k$th leading principle minor of C is given by

[TABLE]

where we have assumed that ${\bf C}_{1},{\bf C}_{2},\cdots,{\bf C}_{n}$ are nonsingular. Using induction, it is not difficult to show that there is a LU decomposition for the correlation matrix. Using the symmetry of C, we write

[TABLE]

which implies that

[TABLE]

The l.h.s of the equation is upper triangular, whereas the r.h.s. is a lower triangular matrix. Consequently, there is a diagonal matrix D such that ${\bf D}={\bf U}({\bf L}^{\rm T})^{-1}$. Then, ${\bf U}={\bf DL}^{\rm T}$, which for the correlation matrix implies that, ${\bf C}={\bf LDL}^{\rm T}$, where D is a positive-definite matrix with its elements also being positive. Accordingly, we write C as ${\bf C}=\tilde{{\bf L}}\tilde{{\bf L}}^{\rm T}$, where $\tilde{{\bf L}}={\bf LD}^{1/2}$, which is the Cholesky decomposition.

It is clear that the matrix $\tilde{{\bf L}}$ is a lower triangular matrix as well, and can be used to transform independent normal variables into dependent multinormal variables, which is the main idea of the method we propose to construct the exact trajectories. The matrix $\tilde{{\bf L}}$ is calculated by [40,41]

[TABLE]

where, $R_{n}=c_{ij}|_{i,j=1}^{n}$ is a positive-definite correlation matrix, $R^{-1}$ is its inverse, and $c_{i}^{*j}=(c_{1j},c_{2j},\dots,c_{i-1j})$ for $j\geq i$, so that $c_{i}\equiv c_{i}^{*i}$. We note that for a semi-positive definite matrix we should remove the first row and first column of the matrix in order to have a positive-definite matrix, and then apply the Cholesky decomposition.

Algorithmically, our Cholesky decomposition algorithm constructs L as follows:

input $n,C_{ij}$

for $k=1,2,...,n$ do

\qquad L_{kk}\leftarrow\Big{(}C_{kk}-\displaystyle\sum_{s=1}^{k-1}L_{ks}^{2}\Big{)}^{1/2}

for $i=k+1,k+2,...,n$ do

\qquad\qquad L_{ik}\leftarrow\Big{(}C_{ik}-\displaystyle\sum_{s=1}^{k-1}L_{is}L_{ks}\Big{)}\Big{/}L_{kk}

end

end

output $L_{ij}$

All the trajectories for FBM in this paper were constructed using this algorithm.

†tabar@uni-oldenburg.de

‡moe@usc.edu

• [1]

S. Redner, A Guide to First-Passage Processes (Cambridge University Press, Cambridge, 2001).

• [2]

S. Condamin, O. Bénichou, V. Tejedor, R. Voituriez, and J. Klafter, First-passage times in complex scale-invariant media, Nature 450, 77 (2007).

• [3]

T. Guérin, N. Levernier, O. Bénichou, and R. Voituriez, Mean first-passage times of non-Markovian random walkers in confinement, Nature 534, 356 (2016).

• [4]

S.A. Rice, Diffusion-Limited Reactions, vol. 25 (Elsevier, Amsterdam, 1985).

• [5]

M.J. Saxton, A biological interpretation of transient anomalous subdiffusion. II. Reaction kinetics, Biophys. J. 94, 760 (2008).

• [6]

B.A. Carreras, V.E. Lynch, I. Dobson, and D.E. Newman, Critical points and transitions in an electric power transmission model for cascading failure blackouts, Chaos 12, 985 (2002).

• [7]

A.L. Lloyd and R.M. May, Epidemiology - how viruses spread among computers and people, Science 292, 1316 (2001).

• [8]

O. Benichou, C. Loverdo, M. Moreau, and R. Voituriez, Two-dimensional intermittent search processes: An alternative to Lévy flight strategies, Phys. Rev. E 74, 020102 (2006).

• [9]

G.M. Viswanathan, E.P. Raposo, and M.G.E. da Luz, Lévy flights and superdiffusion in the context of biological encounters and random searches, Phys. Life Rev. 5, 133 (2008).

• [10]

D. Ben-Avraham and S. Havlin, Diffusion and Reactions in Fractals and Disordered Systems (Cambridge University Press, London, 2000).

• [11]

M. Sahimi, H.T. Davis, and L.E. Scriven, Dispersion in disordered porous media, Chem. Eng. Commun. 23, 329 (1983).

• [12]

M. Sahimi, B.D. Hughes, L.E. Scriven, and H.T. Davis, Dispersion in flow through porous media: I. One-phase flow, Chem. Eng. Sci. 41, 2103 (1986).

• [13]

H.C. Tuckwell, Introduction to Theoretical Neurobiology (Cambridge University Press, London, 1988).

• [14]

A.N. Burkitt, A review of the integrate-and-fire neuron model: I. Homogeneous synaptic input, Biol. Cybern. 95, 1 (2006).

• [15]

L. Sacerdote and M.T. Giraudo, in, Stochastic Biomathematical Models, edited by M. Bachar, J. Batzel, and S. Ditlevsen (Springer, New York, 2013), p. 99.

• [16]

T. Verechtchaguina, I.M. Sokolov, and L. Schimansky-Geier, First passage time densities in non-Markovian models with subthreshold oscillations, Europhys. Lett. 73, 691 (2006).

• [17]

R. Pastor-Satorras and A. Vespignani, Epidemic spreading in scale-free networks, Phys. Rev. Lett. 86, 3200 (2001).

• [18]

O. Benichou, M. Coppey, M. Moreau, P.-H. Suet, and R. Voituriez, Optimal search strategies for hidden targets, Phys. Rev. Lett. 94, 198101 (2005).

• [19]

Q. Hu, Y. Wang, and X. Yang, The hitting time density for a reflected Brownian motion, Comput. Econ. 40, 1 (2012).

• [20]

J. Janssen, O. Manca, and R. Manca, Applied Diffusion Processes, from Engineering to Finance (Wiley, New York, 2013).

• [21]

V. Linetsky, Lookback options and diffusion hitting times: A spectral expansion approach, Finance Stoch. 8, 373398 (2004).

• [22]

D.J. Navarro and I.G. Fuss, Fast and accurate calculations for first-passage times in Wiener diffusion models, J. Math. Psych. 53, 222 (2009).

• [23]

M. Musso and R.K. Sheth, The importance of stepping up in the excursion set approach, Mon. Notices R. Astro. Soc. 438, 2683 (2014).

• [24]

M. Musso and R.K. Sheth, The excursion set approach in non-Gaussian random fields, Mon. Notices R. Astr. Soc. 439, 3051 (2014).

• [25]

V. Pieper, M. Domin, and P. Kurth, Level crossing problems and drift reliability, Math. Methods Oper. Res. 45, 347 (1997).

• [26]

C. Zucca and P. Tavella, The clock model and its relationship with the Allan and related variances, IEEE Trans. Ultras. Ferroelectrics and Frequency Control 52, 289 (2005).

• [27]

C. Zucca, P. Tavella, and G. Peskir, Detecting atomic clock frequency trends using an optimal stopping method, Metrologia 53, S89 (2016).

• [28]

I.M. Sokolov, Models of anomalous diffusion in crowded environments, Soft Matter 8, 9043 (2012).

• [29]

R. Friedrich, J. Peinke, M. Sahimi, and M.R. Rahimi Tabar, Approaching complexity by stochastic methods: from biological systems to turbulence, Phys. Rep. 506 87 (2011).

• [30]

M. Anvari, M.R. Rahimi Tabar, J. Peinke, and K. Lehnertz, Disentangling the stochastic behavior of complex time series, Sci. Rep. 6, 35435 (2016).

• [31]

Q.-H. Wei, C. Bechinger, and P. Leiderer, Single-file diffusion of colloids in one-dimensional channels, Science 287, 625 (2000).

• [32]

T. Turiv, I. Lazo, A. Brodin1, B.I. Lev, V. Reiffenrath, V.G. Nazarenko, and O.D. Lavrentovich, Effect of collective molecular reorientations on Brownian motion of colloids in nematic liquid crystals, Science 342, 1351 (2013).

• [33]

D. Ernst, M. Hellmann, J. Köhler, and M. Weiss, Fractional Brownian motion in crowded fluids, Soft Matter 8, 4886 (2012).

• [34]

T.G. Mason and D.A. Weitz, Optical measurements of frequency-dependent linear viscoelastic moduli of complex fluids, Phys. Rev. Lett. 74, 1250 (1995).

• [35]

F. Nikakhtar, M. Ayromlou, S. Baghram, S. Rahvar, M.R. Rahimi Tabar, and R.K. Sheth, The excursion set approach: Stratonovich approximation and Cholesky decomposition, Mon. Notices R. Astro. Soc. 478, 4, 5296 (2018).

• [36]

T. Franosch, M. Grimm, M. Belushkin, F.M. Mor, G. Foffi, L. Forró, and S. Jeney, Resonances arising from hydrodynamic memory in Brownian motion, Nature 478, 85 (2011).

• [37]

M. Scott, Applied Stochastic Processes in Science and Engineering (Springer, Berlin, 2013).

• [38]

G.R. Jafari, M.S. Movahed, S.M. Fazeli, M.R. Rahimi Tabar, and S.F. Masoudi, Level crossing analysis of the stock markets, J. Stat. Mech. 06, P06008 (2006).

• [39]

T. Verechtchaguina, I.M. Sokolov, and L. Schimansky-Geier, First passage time densities in resonate-and-fire models, Phys. Rev. E 73, 031108 (2006).

• [40]

R. Stratonovich, Topics in the Theory of Random Noise, Vol. 2 (Taylor & Francis, London, 1967).

• [41]

P. Hertz, Uber den gegenseitigen durchschnittlichen Abstand von Punkten, die mit bekannter mittlerer Dichte im Raume angeordnet sind, Mathematische Annalen 67, 387 (1909).

• [42]

B.B. Mandelbrot and J.W. van Ness, Fractional Brownian motions, fractional Gaussian noises and applications, SIAM Rev. 10, 422 (1968).

• [43]

M. Ding and W. Yang, Distribution of the first return time in fractional Brownian motion and its application to the study of on-off intermittency, Phys. Rev. E 5, 207 (1995).

• [44]

J. Krug, H. Kallabis, S.N. Majumdar, S. Cornell, A.J. Bray, and C. Sire, Persistence exponents for fluctuating interfaces, Phys. Rev. E 56, 2702 (1997).

• [45]

G. M. Molchan, Maximum of a fractional Brownian motion: probabilities of small values, Commun. Math. Phys. 205, 97 (1999).

• [46]

J.E. Gentle, Numerical Linear Algebra for Applications in Statistics (Springer, Berlin, 1998), p. 93.

• [47]

V. Madar, Direct formulation to Cholesky decomposition of a general nonsingular correlation matrix, Stat. Probab. Lett. 103, 142 (2015).

• [48]

T. Engel, Firing Statistics in Neurons as Non-Markovian First Passage Time Problem, Ph.D. Thesis, Humboldt-Universität zu Berlin (2006).

• [49]

D.R. Cox and V. Isham, Point Processes (Chapman & Hall, London, 1980).

• [50]

P. Hänggi, P. Talkner, and M. Borkovec, Reaction-rate theory: fifty years after Kramers, Rev. Mod. Phys. 62, 251 (1990).

• [51]

C.R. Doering, P.S. Hagan, and C.D. Levermore, Bistability driven by weakly colored Gaussian noise: The Fokker-Planck boundary layer and mean first-passage times, Phys. Rev. Lett. 59, 2129 (1987).

• [52]

P. Hänggi, P. Jung, and P. Talkner, Comment on “Bistability driven by weakly colored Gaussian noise: The Fokker-Planck boundary layer and mean first-passage times,” Phys. Rev. Lett. 60, 2804 (1988).

• [53]

R. Graham and T. Tél, Nonequilibrium potential for coexisting attractors, Phys. Rev. A 33, 1322 (1986).

• [54]

M.I. Dykman, P.V.E. McClintock, V.N. Smelyanski, N.D. Stein, and N.G. Stocks, Optimal paths and the prehistory problem for large fluctuations in noise-driven systems, Phys. Rev. Lett. 68, 2718 (1992).

• [55]

S.M. Soskin, Large fluctuations in multiattractor systems and the generalized Kramers problem, J. Stat. Phys. 97, 609 (1999).

Table I: The exact enumeration method. The $n$th column corresponds to the $n$th term of the sum in Eq. (5).

[TABLE]