Gap statistics close to the quantile of a random walk

TL;DR

This paper derives the exact distribution of gaps between ordered maxima in a large random walk with Laplace jumps, revealing universal fluctuation behavior and non-universal large deviations.

Contribution

It provides an explicit universal distribution for typical gaps and analyzes their moments, extending understanding of order statistics in correlated stochastic processes.

Findings

Universal alpha-dependent gap distribution with inverse cubic tail

Typical fluctuations scale as O(n^{-1/2})

Large deviations are non-universal

Abstract

We consider a random walk of $n$ steps starting at $x_0=0$ with a double exponential (Laplace) jump distribution. We compute exactly the distribution $p_{k,n}(\Delta)$ of the gap $d_{k,n}$ between the $k^{\rm th}$ and $(k+1)^{\rm th}$ maxima in the limit of large $n$ and large $k$, with $\alpha=k/n$ fixed. We show that the typical fluctuations of the gaps, which are of order $O( n^{-1/2})$, are described by a universal $\alpha$-dependent distribution, which we compute explicitly. Interestingly, this distribution has an inverse cubic tail, which implies a non-trivial $n$-dependence of the moments of the gaps. We also argue, based on numerical simulations, that this distribution is universal, i.e. it holds for more general jump distributions (not only the Laplace distribution), which are continuous, symmetric with a well defined second moment. Finally, we also compute the large deviation form of the gap distribution $p_{\alpha n,n}(\Delta)$ for $\Delta=O(1)$, which turns out to be non-universal.

Figures

Equations

$$x_i=x_{i-1}+{\eta }_i\mathrm{for}i=1,\cdots ,n\mathrm{starting}\mathrm{from}{x}_0=0,$$

$$f(\eta )=\frac{e^{-\sqrt{2}\frac{|\eta |}{\sigma }}}{\sqrt{2}\sigma },$$

$$M_{1,n}=x_{\max }\ge M_{2,n}\ge \cdots \ge M_{n,n}\ge M_{n+1,n}=x_{\min }.$$

$$P_{k,n}(x)=P_{n+2-k,n}(-x).$$

$$P_{k=\alpha n,n}(x)$$

$$\mathrm{with}{P}_{\alpha }(z)$$

$$⟨M_{k,n}⟩\approx \sigma \sqrt{n}\mathcal{M}\left(\alpha =\frac{k}{n}\right),\mathcal{M}(\alpha )=\sqrt{\frac{2}{\pi }}\left(\sqrt{1-\alpha }-\sqrt{\alpha }\right).$$

$$\frac{M_{k=\alpha n,n}}{\sigma \sqrt{n}}\underset{n\to \infty }{⟶}{q}_{1-\alpha }=\inf \{x:{\int }_0^1\Theta (x(\tau )-x)d\tau \ge 1-\alpha \},$$

$$d_{k,n}=M_{k,n}-M_{k+1,n}\ge 0,k=1,\cdots ,n,$$

$$\frac{⟨d_{k,n}⟩}{\sigma }=\frac{\Gamma \left(k+\frac{1}{2}\right)}{\sqrt{2\pi }k!}+\frac{\Gamma \left(n-k+\frac{3}{2}\right)}{\sqrt{2\pi }(n-k+1)!}.$$

$$\frac{⟨d_{k,n}⟩}{\sigma }\approx \{\begin{array}{ll}\frac{\Gamma \left(k+\frac{1}{2}\right)}{\sqrt{2\pi }k!}& ,n\to \infty ,k=O(1)\\ & \\ \frac{\mu (\alpha )}{\sqrt{n}}& ,n\to \infty ,\alpha =\frac{k}{n}=O(1),\end{array}$$

$$\mu (\alpha )=\frac{1}{\sqrt{2\pi }}\left(\frac{1}{\sqrt{\alpha }}+\frac{1}{\sqrt{1-\alpha }}\right).$$

$$p_{k,n}(\Delta )\approx \frac{\sqrt{k}}{\sigma }P\left(\frac{\sqrt{k}\Delta }{\sigma }\right)$$

$$P(\delta )=4\left[\sqrt{\frac{2}{\pi }}(1+2{\delta }^2)-\delta (4{\delta }^2+3)e^{2{\delta }^2}\mathrm{erfc}(\sqrt{2}\delta )\right].$$

$$P(\delta )\approx \{\begin{array}{ll}& 4\sqrt{\frac{2}{\pi }},\delta \to 0,\\ & \\ & \frac{3}{\sqrt{8\pi }}\frac{1}{{\delta }^4},\delta \to \infty .\end{array}$$

$$p_{k=\alpha n,n}(\Delta )\approx \frac{\sqrt{n}}{\sigma }{\mathcal{P}}_{\alpha }\left(\frac{\sqrt{n}\Delta }{\sigma }\right),$$

$${\mathcal{P}}_{\alpha }(\delta )={\int }_0^{\infty }{y}^2{e}^{-\delta y}$$

$$+\frac{y{e}^{-\frac{y^2}{8\alpha }}}{4\sqrt{2\pi }{\alpha }^{\frac{3}{2}}}\mathrm{erfc}\left(\frac{y}{2\sqrt{2(1-\alpha )}}\right)]dy.$$

$${\mathcal{P}}_{\alpha }(\delta )\approx \{\begin{array}{ll}4\sqrt{\frac{2}{\pi }}(\sqrt{\alpha }+\sqrt{1-\alpha }-1)& ,\delta \to 0,\\ & \\ \frac{2}{\pi \sqrt{\alpha (1-\alpha )}}\frac{1}{{\delta }^3}& ,\delta \to \infty .\end{array}$$

$$F_{k,n}(x)=Prob.\left[M_{k,n}\le x\right]=Prob.\left[N_x\le k\right],$$

$$F_{k,n}(x)=\{\begin{array}{ll}\sum _{l=0}^{k-1}{q}_{l,n}(x)& ,x\ge 0\\ & \\ \sum _{l=0}^{k-2}{q}_{n-l,n}(-x)& ,x<0.\end{array}$$

$$q_{k,n}(x)={\int }_0^{\infty }d{x}^{\prime }f(x^{\prime }-x)q_{k,n-1}(x^{\prime })+{\int }_{-\infty }^{0}d{x}^{\prime }f(x^{\prime }-x)q_{n-k,n-1}(-x^{\prime }),$$

$$q_{k,n}(x)={\int }_0^{\infty }d{x}^{\prime }f(x^{\prime }-x)q_{k,n-1}(x^{\prime })+{\int }_0^{\infty }d{x}^{\prime }f(x+x^{\prime })r_{k-1,n-1}(x^{\prime })$$

$$r_{k,n}(x)={\int }_0^{\infty }d{x}^{\prime }f(x^{\prime }-x)r_{k-1,n-1}(x^{\prime })+{\int }_0^{\infty }d{x}^{\prime }f(x+x^{\prime })q_{k,n-1}(x^{\prime })$$

$$\tilde{q}(z,s;x)=\sum _{n=0}^{\infty }\sum _{k=0}^n{s}^n{z}^k{q}_{k,n}(x)\mathrm{and}\tilde{r}(z,s;x)=\sum _{n=0}^{\infty }\sum _{k=0}^n{s}^n{z}^k{q}_{n-k,n}(x),$$

$$\tilde{q}(z,s;x)$$

$$\tilde{r}(z,s;x)$$

$$\frac{{\sigma }^2}{2}{\partial }_x^2\tilde{q}(z,s;x)$$

$$\frac{{\sigma }^2}{2}{\partial }_x^2\tilde{r}(z,s;x)$$

$$\tilde{q}(z,s;x)=a(z,s)e^{-\sqrt{2(1-s)}\frac{x}{\sigma }}+\frac{1}{1-s},\tilde{r}(z,s;x)=b(z,s)e^{-\sqrt{2(1-zs)}\frac{x}{\sigma }}+\frac{1}{1-zs}.$$

$$\tilde{q}(z,s;x)=\left(\frac{1}{\sqrt{1-zs}}-\frac{1}{\sqrt{1-s}}\right)\frac{e^{-\sqrt{2(1-s)}\frac{x}{\sigma }}}{\sqrt{1-s}}+\frac{1}{1-s},$$

$$\tilde{r}(z,s;x)=\left(\frac{1}{\sqrt{1-s}}-\frac{1}{\sqrt{1-zs}}\right)\frac{e^{-\sqrt{2(1-zs)}\frac{x}{\sigma }}}{\sqrt{1-zs}}+\frac{1}{1-zs}.$$

$$\tilde{P}(z,s;x)=\sum _{n=0}^{\infty }\sum _{l=0}^n{s}^n{z}^k{P}_{k,n}(x)=\{\begin{array}{ll}\frac{z}{1-z}{\partial }_x\tilde{q}(z,s;x)& ,x\ge 0\\ & \\ \frac{z^2}{1-z}{\partial }_x\tilde{r}(z,s;-x)& ,x<0.\end{array}$$

Full text

1 Introduction

During the last decades, extreme value statistics (EVS) have found a lot of applications in statistical physics, ranging from disordered systems [1, 2, 3, 4], directed polymers and stochastic growth processes in the Kardar-Parisi-Zhang universality class [5, 6, 7, 8, 9, 10, 11, 12] to fluctuating interfaces [13, 14, 15, 16], random matrices [17, 18], random walks and Brownian motions [19, 20, 21] all the way to cold atoms [22]. The basic question concerns the distribution of the maximum $x_{\max}$ (or equivalently the minium $x_{\min}$) among a collection of $N$ random variables $x_{1},x_{2},\cdots,x_{N}$, and in particular in the limit $N\to\infty$, i.e. in the thermodynamic limit. This problem is fully understood in the case of independent and identically distributed (i.i.d.) random variables $x_{i}$’s, for which it is well known that there exist three distinct universality classes (Gumbel, Fréchet and Weibull) depending only on the tail of the parent distribution of the $x_{i}$’s [23]. However, in many situations in statistical physics, it turns out that, often, one has to deal with strongly correlated variables [24]. In fact, there exist at present very few exact results for the EVS in strongly correlated systems and it is thus crucial to identify physically relevant models for which the EVS can be computed exactly. A prototypical example of such models is the discrete-time random walk (RW), which constitutes a useful laboratory to test the effects of strong correlations on EVS [25]. Here we are interested in the statistics of the gaps between the consecutive maxima of a discrete-time RW.

Indeed, the distribution of the global maximum $x_{\max}$ (or the minimum $x_{\min}$) is certainly interesting but it gives only a partial information on the system – it concerns one single variable out of $n\gg 1$ – and in some cases it is useful to consider the more general question of order statistics which concern the joint statistics of the $k$-th maxima $M_{k,n}$ such that $x_{\max}=M_{1,n}>M_{2,n}>\cdots>M_{n+1,n}=x_{\min}$. Natural observables are then the gaps between successive maxima, $d_{k,n}=M_{k,n}-M_{k+1,n}$, which are useful, for instance, to quantify the phenomenon of “crowding” near the extremes [26, 27, 28, 29]. In physics, the statistics of the gaps were studied in the context of branching Brownian motions [30, 31], as well as for noisy signals with power spectrum in $1/f^{\alpha}$ [32], with applications in cosmology [33]. For RWs, such questions related to the $k$-th maximum belong to the general realm of “fluctuation theory” [34] and their statistics have been computed using probabilistic methods [35, 36, 37, 38, 39, 40, 41]. However, much less is known about the gaps from fluctuation theory (see however [42]).

Rather recently, two of us developed an independent method, based on backward equations (see below), which allowed us to solve exactly the gap statistics for random walks with symmetric exponential jumps [20]. From this exact result, the distribution of the gaps near $x_{\max}$ (i.e., at the “edge”), for long random walks, was obtained and it was shown to exhibit a very rich behaviour, which was conjectured to be, to a large extent, universal, i.e. independent of the details of the jump distribution provided it has a well defined second moment. This conjecture at the “edge” was partly confirmed by a more recent work where it was shown that, for long RWs, the same distribution describes the typical fluctuations of the gaps of RWs with symmetric gamma-distributed jumps [43].

The goal of this paper is to reconsider this problem of the gaps of random walks and study their statistics in the “bulk”, i.e. far away from $x_{\max}$ 111Note that the terms “bulk” and “edge” are borrowed from the terminology used in Random Matrix Theory [44, 45]., near a quantile of the random walk, for instance near the median which is at half-way between $x_{\min}$ and $x_{\max}$. We find that the gaps in the bulk also display a very rich behaviour which is quite different from the one found at the edge. We conjecture that this behaviour is also universal.

2 Model and main results

Let us thus consider a one-dimensional random walk in continuous space and discrete time defined by

[TABLE]

where the $\eta_{i}$’s are i.i.d. random variables with a double exponential (or Laplace) probability distribution function (PDF)

[TABLE]

where $\sigma^{2}=\int_{-\infty}^{\infty}\eta^{2}\,f(\eta)\,d\eta$ is the variance of the jump distribution. Hence in the large $n$ limit the random walk in (1) converges to the Brownian motion. Note that for a walk of $n$ steps, there are $n+1$ positions $\{x_{0},x_{1},\cdots,x_{n}\}$. We order these positions and define the random variables $M_{k,n}$ as the $k^{\rm th}$ maximum among the positions $x_{i}$’s of the random walk (see Fig. 1), such that

[TABLE]

Since the jump distribution is symmetric, i.e. $f(\eta)=f(-\eta)$, $x_{\max}=M_{1,n}$ has the same distribution as $-x_{\min}=-M_{n+1,n}$. Similarly $M_{2,n}$ has the same distribution as $-M_{n,n}$ and more generally $M_{k,n}$ has the same distribution as $-M_{n+2-k,n}$. Therefore the distribution $P_{k,n}(x)$ of $M_{k,n}$ satisfies the relation

[TABLE]

Distribution of the $k^{\rm th}$ maximum. We first compute the PDF $P_{k,n}(x)$ of $M_{k,n}$ for a double exponential jump PDF, using the method introduced in [20]. In the limit of large $n$, with $\alpha=k/n$ fixed, one finds that $P_{k=\alpha n,n}(x)$ takes the scaling form

[TABLE]

Note that the limiting distribution satisfies the relation $P_{\alpha}(z)=P_{1-\alpha}(-z)$ which reflects the symmetry between maxima and minima noted in Eq. (4). In Fig. 2, we plot the PDF $P_{\alpha}(x)$ for $\alpha=0.1$. We see clearly that it is quite asymmetric, and from the exact expression (5), it is easy to check that $P_{\alpha}(z)\approx e^{-z^{2}/(2(1-\alpha))}$ for $z\to+\infty$ while $P_{\alpha}(z)\approx e^{-z^{2}/(2\alpha)}$ for $z\to-\infty$. Close to $z=0$, where $P_{\alpha}(0)=\sqrt{2/\pi}$, the distribution $P_{\alpha}(z)$ has a cusp. From this result (5), one can easily compute the first moment $\langle M_{k,n}\rangle$ in the limit of large $n$ and $k$ with $k/n=\alpha$ fixed

[TABLE]

Note that ${\cal M}(1-\alpha)=-{\cal M}(\alpha)$, in agreement with the property in (4). In the limit $\alpha\to 0$, one recovers that $\langle x_{\max}\rangle=\langle M_{1,n}\rangle\approx\sigma\sqrt{n}\,{\cal M}(0)=\sigma\sqrt{2\,n/\pi}$, as expected from the Brownian motion result.

Incidentally, in this scaling limit where both $n$ and $k$ are large with $k/n=\alpha$ fixed, $M_{k,n}$ corresponds exactly to what is called, in the probability literature, the $(1-\alpha)$-quantile for Brownian motion [37]. Roughly speaking, $M_{k=\alpha n,n}$ is such that a fraction $(1-\alpha)$ of the points of the trajectory of the random walks are below $M_{k=\alpha n,n}$, while a fraction $\alpha$ of the points are above it. Since in the large $n$ limit the RW converges to Brownian motion, $M_{k=\alpha n,n}/(\sigma\sqrt{n})$, converges, in the scaling limit where both $n$ and $k$ are large with $k/n=\alpha$ fixed, to the $(1-\alpha)$-quantile $q_{1-\alpha}$ of the Brownian motion, i.e.

[TABLE]

where $x(\tau)$ is a standard Brownian motion (with diffusion constant $D=1/2$) starting from [math] on the time interval $[0,1]$ and ${\Theta(x)}$ is the Heaviside theta function. In fact, one can check that the formula for $P_{\alpha}(z)$ in Eq. (5), obtained here using a backward-equation formalism, coincides with the result obtained previously in the mathematics literature using quite different probabilistic methods for Brownian motion [37, 38, 39].

Distribution of the $k^{\rm th}$ gap. Our main new results concern the gaps

[TABLE]

between two consecutive maxima. For a double exponential jump distribution (2), the mean value of the gap $\langle d_{k,n}\rangle=\langle M_{k,n}\rangle-\langle M_{k+1,n}\rangle$ can be computed for any finite $n$ and $k$, yielding [20]

[TABLE]

It is straightforward to extract the large $n$ behaviour of this exact expression (9) in the two scaling regimes corresponding to $k={O}(1)$ and $k={O}(n)$ as

[TABLE]

where the scaling function $\mu(\alpha)$ reads

[TABLE]

Note that one can check that $\mu(\alpha)=-{\cal M}^{\prime}(\alpha)$, where ${\cal M}(\alpha)$ is given in Eq. (6), as expected from Eqs. (8) and (6). This result (10) clearly shows that, for large $n$, there are two different scales for the gaps $d_{k,n}$ depending on $k=O(1)$ or $k=O(n)$. It is useful to think about the values of the $k$-th maxima of the random walks after step $n$ as a point process on the line, as illustrated in Fig. 3.

Near the edges, i.e. near the maximum $x_{\max}=M_{1,n}$ and the minimum $x_{\min}=M_{n+1,n}$, the gaps are of order $O(1)$ (see the first line of Eq. (10)) while they are of order $O(n^{-1/2})$ (see the second line of Eq. (10)) in the bulk, i.e. far from $x_{\max}$ and $x_{\min}$. Note that by taking the large $k$ limit in the first line of Eq. (10) one obtains $\langle d_{k,n}\rangle\approx\sigma/\sqrt{2\pi k}$, while by taking the small $\alpha=k/n$ limit in the second line of Eq. (10) one also obtains $\langle d_{k,n}\rangle\approx\sigma/\sqrt{2\pi\alpha n}=\sigma/\sqrt{2\pi k}$: this shows that there is a smooth matching between the edge and the bulk at the level of the first moment $\langle d_{k,n}\rangle$.

What about the full PDF $p_{k,n}(\Delta)$ of the gaps $d_{k,n}$ in the large $n$ limit? Near the edge, this PDF was computed for jumps with a double exponential jump distribution in Ref. [20] and subsequently for symmetric gamma-distributed jumps in Ref. [43]. In particular, it was shown in Ref. [20] that for a double exponential jump distribution, the PDF $p_{k,n}(\Delta)$ becomes independent of $n$ in the large $n$ limit, i.e. $\lim_{n\to\infty}p_{k,n}(\Delta)=p_{k,\infty}(\Delta)$, consistent with the first line of Eq. (10). In the large $k$ limit, it turns out [20] that the limiting distribution $p_{k,\infty}(\Delta)$ has two different scaling behaviours, depending of $\Delta$: (i) a regime of typical fluctuations for $\Delta=O(1/\sqrt{k})$ and (ii) a large deviation regime for $\Delta=O(1)$. The most interesting result obtained in [20] concerns the typical fluctuations where $p_{k,n}(\Delta)$ takes the scaling form [20]

[TABLE]

where the scaling function $P(\delta)$ is given by [20]

[TABLE]

Based on numerical simulations, it was conjectured in [20] that the typical distribution $P(\delta)$ is universal, i.e. it does not depend on the jump distribution $f(\eta)$ as long as it is symmetric and has a finite variance $\sigma^{2}<\infty$. The validity of this conjecture was then reinforced by an exact analytical computation for gamma distributed jump distribution $f(\eta)=\frac{|\eta|^{p}}{2p!}e^{-|\eta|}$ with $p\in\mathbb{N}$ [43]. From this expression (13), it is easy to obtain the asymptotic behaviours of $P(\delta)$ for small and large $\delta$

[TABLE]

In particular, it exhibits an interesting power law tail $P(\delta)\propto\delta^{-4}$ for large $\delta$. The large deviation regime of $p_{k,n}(\Delta)$, for $\Delta=O(1)$, can also be computed explicitly for the double exponential jump distribution (2) but, unlike $P(\delta)$ in (14), it turns out to be non-universal, i.e. it depends explicitly on the jump distribution [20, 43] and, for this reason, it is somewhat less interesting than the typical fluctuation regime.

In this paper, we derive the full PDF $p_{k,n}(\Delta)$ of the gap $d_{k,n}$ in the bulk, i.e. for large $n$ and large $k$ but keeping the ratio $\alpha=k/n$ fixed. We show that the behaviour in the bulk is rather different from the one found at the edge in [20, 43] recalled above in Eqs. (13) and (14), which corresponds instead to the limit $\alpha\to 0$ (i.e. $1\ll k\ll n$). We find that this PDF $p_{k,n}(\Delta)$ again exhibits two different scaling regimes depending on $\Delta$: a typical regime for $\Delta=O(n^{-1/2})$, consistent with the second line of Eq. (10), and a large deviation regime for $\Delta=O(1)$. Our most interesting results concern the typical regime, for $\Delta=O(n^{-1/2})$, where $p_{k,n}(\Delta)$ takes the scaling form

[TABLE]

where the scaling function ${\cal P}_{\alpha}(\delta)$ depends continuously on the parameter $\alpha$ and is given explicitly by

[TABLE]

For generic $\alpha$ this integral over $y$ can be evaluated in terms of hypergeometric functions of two variables (namely Humbert series, see Eq. (89) in Appendix A). For the special case $\alpha=1/2$, which describes the gaps near the median, ${\cal P}_{1/2}(\delta)$ can be expressed in terms of elementary functions [see Eq. (A)]. One can also show that in the limit $\alpha\to 0$, our result in Eqs. (15) and (16) yields back the edge result in Eqs. (13) and (14) – this limit is however a bit subtle and is studied in detail below. Under this form (16), we notice the relation ${\cal P}_{\alpha}(\delta)={\cal P}_{1-\alpha}(\delta)$, which is a direct consequence of the symmetry observed above for the distribution of the $k$-th maximum in (4). In Fig. 4, we show a plot of this scaling function ${\cal P}_{\alpha}(\delta)$, for $\alpha=1/2$, and compare it to numerical data (appropriately scaled according to Eq. (15)) obtained by simulating random walks (1) for three different jump distributions: a double exponential PDF (2) – for which our computation is exact – but also a Gaussian PDF as well as a uniform PDF, which we can not study analytically. We find that the agreement between our theoretical results (16) and the numerical simulations are equally good for the three jump distributions. Based on this observation, we conjecture that this distribution ${\cal P}_{\alpha}(\delta)$ is universal, i.e. it is independent of the jump distribution $f(\eta)$ provided it is continuous, symmetric and with a well defined variance $\sigma^{2}$.

From this integral representation (16) one can rather easily extract the asymptotic behaviours of the scaling function ${\cal P}_{\alpha}(\delta)$

[TABLE]

Interestingly, we see that the tail ${\cal P}_{\alpha}(\delta)\propto\delta^{-3}$, for finite $\alpha$ and in the bulk, is different from the tail $\propto\delta^{-4}$ obtained at the edge [see Eq. (14)]. From this inverse cubic tail one would naively conclude that the moments of the gaps (beyond the first one given in Eqs. (10) and (11)) is not defined. However, this power law behaviour of the gap distribution $p_{k,n}(\Delta)$ is cut-off for $\Delta\gg n^{-1/2}$ and the higher moments are actually dominated by the large deviation regime of the PDF $p_{k,n}(\Delta)$ for $\Delta=O(1)\gg n^{-1/2}$ which we can also compute exactly (see below). The latter turns out to be non-universal. Consequently, the moments of the gaps beyond the first one, that we also study below, are to a large extent non-universal.

The paper is organised as follows. Section 3 is dedicated to the distribution of the $k^{\rm th}$ maximum $M_{k,n}$. In Section 4 we derive the results for the gap $d_{k,n}$, which constitute our main results, before we conclude in Section 5. Some details of the computations have been relegated in Appendices A, B and C.

3 Distribution of the $k^{\rm th}$ maximum

We first expose a method which allows us to obtain the distribution of $M_{k,n}$ that we will generalise in the next section to obtain the distribution of the gaps $d_{k,n}$. It relies on the identity for the cumulative distribution function (CDF) of the $k^{\rm th}$ maximum $M_{k,n}$

[TABLE]

where $N_{x}$ is the counting process for the number of steps where the random walk takes values above $x$. Indeed, there are exactly $k$ positions among the $n+1$ positions $x_{i}$’s of the walk such that $x_{i}\leq M_{k,n}$ (see Fig. 1 for an example). Note that this identity remains valid for any discrete time stochastic process. We introduce the probability $q_{k,n}(x)$ that a random walk of $n$ steps starting at $x_{0}=x$ has exactly $k$ points on the negative axis between step $1$ and step $n$, i.e. $N_{0}=n-k$ for the walk starting from $x_{0}=x$. The CDF $F_{k,n}(x)$ of $M_{k,n}$ can then be expressed from an elementary path transformation (see Fig. 5) as [20]

[TABLE]

The probability $q_{k,n}(x)$ can be constructed recursively, using the equation [20]

[TABLE]

together with the initial condition $q_{0,0}(x)=1$ and $q_{k,n}(x)=0$ for $k>n$. The first term of this equation describes the case where the walk has an additional initial jump from $x>0$ to $x^{\prime}>0$, while the second term describes a jump from $x>0$ to $x^{\prime}<0$. To solve this equation (20) it is useful [20] to introduce the auxiliary function $r_{k,n}(x)=q_{n-k,n}(x)$ which is the probability that a RW starting from $x_{0}=x$ has $k$ points above [math] between step $1$ and step $n$. One can then write two coupled equations for $q_{k,n}(x)$ and $r_{k,n}(x)$ [20]

[TABLE]

These integral equations (21) can be solved by generating function techniques. We introduce

[TABLE]

and obtain from (21) the set of coupled integral equations

[TABLE]

These equations, valid for any distribution of jumps $f(\eta)$, turn out to be very difficult to solve in general. However, for a double exponential jump distribution (2) they can be solved exactly using the identity $f^{\prime\prime}(\eta)=\frac{2}{\sigma^{2}}\left[f(\eta)-\delta(\eta)\right]$. Differentiating twice Eqs. (24) and (25) with respect to $x$, we obtain two decoupled differential equations

[TABLE]

Discarding the diverging solution for $x\to+\infty$, we obtain

[TABLE]

The values of $a(z,s)$ and $b(z,s)$ are obtained by substituting back these forms in Eqs. (26) and (27). Finally, the solutions read

[TABLE]

The generating function $\tilde{P}(z,s;x)$ of the PDF $P_{k,n}(x)=F_{k,n}^{\prime}(x)$ can be worked out explicitly in terms of $\tilde{q}(z,s;x)$ and $\tilde{r}(z,s;x)$ using Eq. (19)

[TABLE]

To obtain the large $n$ behaviour of $P_{k,n}(x)$, we perform a change of variables in the generating function $\tilde{P}(z,s;x)$ and set $s=\exp(-p)$ and $z=\exp(-q)$ with $p\sim q\ll 1$. In this limit, the discrete sums over $n$ and $k$ can be replaced by integrals. The generating function $\tilde{P}(z,s;x)$ therefore converges towards the double Laplace transform of the PDF $P_{k,n}(x)$ with respect to $n$ and $k$,

[TABLE]

where the function $\tilde{\Pi}(q,p;x)$ can be computed from Eqs. (29) and (30). It reads, at leading order for $p\sim q\ll 1$,

[TABLE]

Using the inverse Laplace transforms

[TABLE]

we invert the Laplace transform from $p$ to $n$, yielding

[TABLE]

Finally, using the identity

[TABLE]

we invert the Laplace transform from $q$ to $k$, yielding

[TABLE]

Changing the variable $\tau\to n\tau$ in this expression, it takes the scaling form described in the first line of Eq. (5), where the scaling function $P_{\alpha}(z)$ is given as

[TABLE]

The integral in Eq. (37) can be computed exactly using the identity

[TABLE]

where $\operatorname{erfc}(x)=\frac{2}{\sqrt{\pi}}\int_{x}^{\infty}e^{-u^{2}}du$ is the complementary error function, leading to the final expression of $P_{\alpha}(z)$ in the second line of Eq. (5).

We conclude this section by mentioning that there is an alternative method to obtain the PDF of $M_{k,n}$ for a discrete time random walk. It was obtained in the limit of large $n$ and for $\alpha=k/n=O(1)$ in Ref. [40], making use of the identity in law derived in [35] and extended in [36] (see also the more recent work [41])

[TABLE]

where $\{x_{i}\}$ and $\{x_{i}^{\prime}\}$ are two independent random walks with same jump distribution $f(\eta)$ starting at $x_{0}=x_{0}^{\prime}=0$. This alternative method can be exploited for any jump distribution $f(\eta)$, even if $\sigma^{2}=\infty$. The distribution of $M_{k,n}$ was shown to take universal scaling forms, depending on the behaviour of $\hat{f}(k)=\int_{-\infty}^{\infty}\frac{dk}{2\pi}e^{-ik\eta}f(\eta)$ for small $k$. However, to our knowledge, there is no direct extension of this method to compute the distribution of the gap $d_{k,n}$. In the next section, we will show how to extend the method presented in this section to obtain exact results for the PDF $p_{k,n}(\Delta)$ of the gap $d_{k,n}$ for a double exponential jump distribution.

4 Distribution of the $k^{\rm th}$ gap

Our starting point is the joint CDF of $M_{k,n}$ and $M_{k+1,n}$,

[TABLE]

from which the PDF of the gap $d_{k,n}$ can be obtained as

[TABLE]

To obtain the joint CDF $S_{k,n}(x,y)$, we introduce the probability $Q_{k,n}(x,\Delta)$ that a random walk of $n$ steps starting at $x_{0}=x$ has exactly $k$ points below [math] and no point in the interval $[-\Delta,0]$ between step $1$ and step $n$. Using a simple path transformation, one obtains the relation [20]

[TABLE]

The probability $Q_{k,n}(x,\Delta)$ can be obtain recursively using the relation

[TABLE]

This is a similar recursion relation as for $q_{k,n}(x)$ in Eq. (20). Therefore, introducing the generating functions

[TABLE]

they satisfy a set of coupled integral equations very similar to Eqs. (24) and (25) [20]

[TABLE]

For a double exponential jump distribution $f(\eta)$ as in Eq. (2), these integral equations can be recast as differential equations [20, 43]: indeed one can easily show that the generating functions follow the set of differential equations (26) and (27), with the substitutions $\tilde{q}\to\tilde{Q}$ and $\tilde{r}\to\tilde{R}$. Solving these equations, we obtain

[TABLE]

The coefficients $A_{1}$ and $B_{1}$ are then determined by inserting these solutions back in Eqs. (45) and (46). This yields

[TABLE]

and $B_{1}(z,s;\Delta)=A_{1}(z^{-1},zs;\Delta)$. From Eqs. (41) and (42), we can express the generating function $\tilde{p}(z,s;\Delta)=\sum_{n=0}^{\infty}\sum_{k=0}^{n}s^{n}z^{k}p_{k,n}(\Delta)$, in terms of the coefficients $A_{1}$ and $B_{1}$ (see Appendix A of Ref. [43] for more details). This yields

[TABLE]

As in the case of the PDF of $M_{k,n}$, we are interested in the limit $n\to\infty$ and $k\to\infty$, which is conveniently obtained by performing the changes of variables $z=e^{-q}$ and $s=e^{-p}$ and by taking the limit $p,q\to 0$. In this limit, the discrete sums over $k$ and $n$ can then be replaced by integrals, yielding

[TABLE]

where $\tilde{\pi}(p+q,p;\Delta)$ is the double Laplace transform of the PDF $p_{k,n}(\Delta)$ with respect to $k$ and $n-k$. To simplify the notations, we denote from now on $r=p+q$. Note that we anticipate that $\tilde{\pi}(r=p+q,p;\Delta)$ is symmetric in $p$ and $p+q$, since $p_{k,n}(\Delta)=p_{n-k,n}(\Delta)$. We will now analyse the PDF in the large $n$ limit and treat separately the typical fluctuations for $\Delta=O(n^{-1/2})$ and the atypically large fluctuations for $\Delta=O(1)$.

4.1 Typical regime of fluctuation

To analyse the typical regime, we need to obtain the behaviour of $A_{1}(z=e^{-q},s=e^{-p};\Delta)$ (resp. $B_{1}$) in the regime $p\sim q\sim\Delta^{2}$. In this regime, the coefficients take the scaling form

[TABLE]

Note that $b_{1}(r,p;\Delta)=a_{1}(p,r;\Delta)$. Inserting these expressions in Eq. (50), we realise that the leading terms are the second derivatives with respect to $\Delta$, as $\Delta$ is small and $\sqrt{1-s}\sim\sqrt{p}\ll 1$. The scaling function $\tilde{\pi}(r=p+q,p;\Delta)$ thus reads in this limit

[TABLE]

It is symmetric in $p$ and $r=p+q$, reflecting the symmetry of the PDF $p_{k,n}(\Delta)$ in $k$ and $n-k$. To invert the Laplace transforms with respect to $r$ and $s$ we first use the identity

[TABLE]

to obtain

[TABLE]

Finally, using the Laplace inversion formulae

[TABLE]

we obtain the PDF

[TABLE]

Performing the change of variable $y\to y/\sqrt{n}$ in Eq. (60), we eventually obtain that $p_{k,n}(\Delta)$ takes the scaling form announced in Eq. (15) with the scaling function given in Eq. (16) that we reproduce here

[TABLE]

For generic $\alpha$ this integral has a rather complicated expression in terms of hypergeometric functions of two variables (namely Humbert series, see Eq. (89) in Appendix A). However, for $\alpha=1/2$ it has an explicit expression in terms of elementary functions given in (A).

From this expression (4.1), we can compute the mean of the distribution, recovering the result of Eq. (11) (see Appendix B for details). However, this distribution has a heavy tail as seen in Eq. (17) and its moments of order $p\geq 2$ are infinite. In Fig. 4, we compare the scaling function ${\cal P}_{\alpha}(\delta)$ to numerical results obtained for $10^{6}$ simulations of random walks of $n=10^{3}$ steps with exponential, Gaussian and uniform distribution of jump $f(\eta)$, suggesting the universality of the result.

The limit $\alpha\to 0$. We first check that in the limit $\alpha\to 0$, this distribution ${\cal P}_{\alpha}(\delta)$ yields back the result at the edge obtained in Ref. [20], given in Eqs. (12) and (13). To recover this edge result, we need to take simultaneously the limit $\alpha=k/n\to 0$ and $\delta=\sqrt{n}\Delta\to\infty$ but keeping $\sqrt{\alpha}\,\delta=\sqrt{k}\,\Delta$ fixed [see Eq. (13)]. In this scaling limit, we show that ${\cal P}_{\alpha}(\delta)$ in Eq. (16) takes the scaling form

[TABLE]

where $P(\delta)$ is given in Eq. (13). To show this result (62), we demonstrate equivalently, setting $\delta=\delta^{\prime}/\sqrt{\alpha}$ with $\delta^{\prime}$ fixed, that

[TABLE]

To show (63) we write ${\cal P}_{\alpha}(\delta^{\prime}/\sqrt{\alpha})$ starting from Eq. (4.1) and perform the change of variable $y\to z=y/\sqrt{\alpha}$ to obtain

[TABLE]

since only the last term (in the integrand) survives in the limit $\alpha\to 0$. Finally, evaluating explicitly the remaining integral over $z$ in Eq. (64) yields back the scaling function $P(\delta^{\prime})$ given in Eq. (13). This shows the scaling form in Eq. (62).

The limit $\alpha\to 0$ deserves yet another remark. Indeed, in the limit $\delta\to\infty$, the scaling function ${\cal P}_{\alpha}(\delta)$, in the bulk, has an inverse cubic tail [see the second line of Eq. (17)], i.e. ${\cal P}_{\alpha}(\delta)\propto\delta^{-3}$. On the other hand, at the edge (corresponding to the limit $\alpha\to 0$), the PDF of the gap decays as $P(\delta)\propto\delta^{-4}$ [see Eq. (14)]. This indicates that the two limits $\alpha\to 0$ and $\delta\to\infty$ do not commute. In fact, from the full expression in (4.1) it is rather straightforward to show that there exists a scaling regime corresponding to $\alpha\to 0$ and $\delta\to\infty$ but keeping $\xi=\alpha\delta=O(1)$ fixed, which smoothly interpolates between these two different tail behaviours. In this scaling regime, we find that ${\cal P}_{\alpha}(\delta)$ takes the scaling form

[TABLE]

A plot of this scaling function ${\cal F}(\xi)$ together with a comparison of $\alpha^{-5/2}{\cal P}_{\alpha}(\delta)$ [evaluated from the exact formula in Eq. (4.1)] is provided in the left panel of Fig. 6. This scaling form (65) indicates that, for large $\delta\gg 1$ and small $\alpha\ll 1$, the function ${\cal P}_{\alpha}(\delta)$ exhibits two different tail behaviours: if $1\ll\delta\ll\alpha^{-1}$ then ${\cal P}_{\alpha}(\delta)\sim\delta^{-4}$ while if $\delta\gg\alpha^{-1}$ then $P_{\alpha}(\delta)\propto\delta^{-3}$. This is summarised in the right panel of Fig. 6.

4.2 Large deviation regime

Since the PDF governing the typical fluctuations of the gaps has a power law tail $P_{\alpha}(\delta)\propto\delta^{-3}$, higher order moments of the gaps are dominated by the large deviation regime of $p_{k,n}(\Delta)$ for $\Delta=O(1)$. To study this regime, we compute the behaviour of the coefficients $A_{1}(z=e^{-q},s=e^{-p};\Delta)$ for $p\sim q\ll 1$ but with $\Delta=O(1)$ fixed. This yields

[TABLE]

where the functions $a_{2}(r=p+q,p;\tilde{\Delta})$ and $a_{3}(r=p+q,p;\tilde{\Delta})$ read

[TABLE]

In the limit of small $p\sim q\ll 1$, the leading contribution in the general formula (50) is again given by the terms involving the second derivatives, yielding

[TABLE]

In this expression, we need to keep only the terms that depend both on $r$ and $p$ as they are the only terms giving physical contribution when inverting the Laplace transforms. These terms read

[TABLE]

For this large deviation form, the Laplace transforms are simple to invert, using that

[TABLE]

For large $n$ and large $k$ we obtain finally for $\Delta=O(1)$

[TABLE]

where $\Psi(\tilde{\Delta})$ and $\varphi_{0}(\tilde{\Delta})$ are given by

[TABLE]

The leading behaviour of this large deviation form (73) will be different in the large $n$ limit, depending on whether $k=O(1)$ remains fixed or $k=O(n)$. Indeed, taking first the limit $n\to\infty$ with $k$ fixed (i.e. $\alpha\to 0$, corresponding to the edge), the dominant contribution is the term $\propto\alpha^{-3/2}$ in Eq. (74), leading to

[TABLE]

recovering the result of Ref. [20].

On the other hand, in the bulk, with $n,k\to\infty$ with $\alpha=k/n$ fixed, the leading contribution is given by the term of order $O(1/n)$ in (73). Therefore, in the bulk, the large deviation form of the PDF of the gap reads

[TABLE]

Note that in the limit $\tilde{\Delta}\to 0$, the large deviation function $\Psi(\tilde{\Delta})$ behaves as $\Psi(\tilde{\Delta})\approx(2/\pi)\tilde{\Delta}^{-3}$, which matches smoothly with the tail behaviour of the typical regime [see the second line of Eq. (17)].

Finally, we end up this section on the large deviations by noting that all the terms in Eq. (73) become of the same order in the intermediate regime where $k=O({\sqrt{n}})$. Indeed, setting $\lambda=k/\sqrt{n}=O(1)$ one obtains, from (73), that $p_{k,n}(\Delta)$ takes the scaling form

[TABLE]

which smoothly interpolates between (75) in the limit $\lambda\to 0$ and (76) in the limit $\lambda\to\infty$.

Let us now investigate the consequences of this behaviour (73) on the moments of the gaps.

4.3 Computation of the moments

Since the scaling function ${\cal P}_{\alpha}(\delta)$ that describes the typical gap fluctuations behaves as ${\cal P}_{\alpha}(\delta)\sim\delta^{-3}$ for large $\delta$, the first moment of the gap is indeed completely dominated by the typical region [see Eq. (96)]. This is however not the case for higher order moments. Indeed, the two non-trivial contributions to the large deviation form in Eq. (73) will both contribute to the moments of order $p\geq 2$. In the case $p>3$, we obtain

[TABLE]

where the values of $M_{p}$ and $m_{p}$ can be computed explicitly (see Appendix C),

[TABLE]

where $\zeta(s)=\sum_{k=1}^{\infty}k^{-s}$ is the Riemann Zeta function. The cases $p=2$ and $p=3$ are particular since there are logarithmic corrections. For $p=2$ we obtain

[TABLE]

while, for $p=3$ we have

[TABLE]

Finally, these behaviour can be summarised as follows

[TABLE]

In the regime $\alpha=k/n=O(1)$, the first term in the last three lines of Eq. (83) gives the leading contribution to the moments. On the other hand, in the regime $k=O(1)$, it is the term in $k^{-3/2}$ that is dominant. Finally, in the intermediate regime mentioned above $k=O(\sqrt{n})$, both of these terms are of the same order for $p>3$ while the term with the logarithmic correction is dominant for $p=2,3$.

5 Conclusion

In this article, we have computed exactly the PDF $p_{k,n}(\Delta)$ of the gap between two successive maxima $M_{k,n}$ and $M_{k+1,n}$ for a random walk with double exponential (Laplace) jump distribution. The main focus of the present paper has been the limiting distribution of the gap $d_{k,n}$ in the scaling limit where both $n$ and $k$ are large, keeping the ratio $\alpha=k/n$ fixed. This allowed us to study the gaps in the bulk, i.e. far from the global maximum $x_{\max}$ of the random walk after $n$ steps (see Fig. 3). Our main result is an explicit expression for the distribution ${\cal P}_{\alpha}(\delta)$ [see Eq. (16)] which governs the typical fluctuations of $d_{k,n}$ in this scaling limit, namely for $d_{k,n}=O(n^{-1/2})$. We conjecture that this distribution ${\cal P}_{\alpha}(\delta)$ is universal for all random walks with a jump distribution $f(\eta)$ which is continuous, symmetric and possesses a finite second moment. What happens for heavy tailed jump distributions, i.e. the case of Lévy flights, remains a challenging open question, in particular because we do not know how to solve the backward integral equations (45) and (46) in this case. We hope that the results obtained here will motivate further works to develop alternative methods to study the gap statistics of Lévy flights.

We found rather useful to think about the different positions of the random walker after $n$ steps as a point process on a line, as illustrated in Fig. 3. By analogy with random matrices we naturally identify edge regions, close to the extremal positions of the random walk, as well as the bulk region, far from the maximum and the minimum. In particular, we have shown that the gaps behave quite differently in these two regions. Pursuing this analogy with random matrices, one may wonder whether one can define a “density” associated to this point process that would capture the existence of these edge regions, and would be the equivalent of the Wigner semi-circle in random matrix theory. This is left for future investigations [46].

Appendix A Explicit expression for ${\cal P}_{\alpha}(\delta)$

The integrals in the expression (16) for the PDF ${\cal P}_{\alpha}(\delta)$ can be expressed in terms of the following integrals

[TABLE]

as well as (see formula 4 p. 178 of Ref. [47])

[TABLE]

where $\Psi_{1}$ is a confluent hypergeometric series of two variables (sometimes called Humbert series [48]) defined as

[TABLE]

where $(a)_{m}=\Gamma(a+m)/\Gamma(a)$ is the Pochhammer symbol. In terms of $I_{2}(\delta,b)$ and $I_{3}(\delta,b)$ in (84) and $J(\delta,b,c)$ in (86), the gap distribution ${\cal P}_{\alpha}(\delta)$ in (16) reads

[TABLE]

Note that this expression (89) is explicitly symmetric under the change $\alpha\to 1-\alpha$, as it should, i.e. ${\cal P}_{\alpha}(\delta)={\cal P}_{1-\alpha}(\delta)$.

In the special case $\alpha=1/2$ (which corresponds to the vicinity of the median), the integrals in Eq. (89) can be performed in terms of elementary functions (using in particular formula 2 p. 175 of [47]). This yields

[TABLE]

which is clearly different from the scaling function found at the edge [20] (corresponding to the limit $\alpha\to 0$), given in Eq. (13). In particular, its asymptotic behaviours are given by

[TABLE]

which is fully consistent with the behaviours given in Eq. (17) specified for $\alpha=1/2$.

Appendix B Computation of $\langle d_{k,n}\rangle$ from ${\cal P}_{\alpha}(\delta)$

From Eq. (4.1), the mean value of the gap is obtained as

[TABLE]

Using that $\int_{0}^{\infty}y^{2}\delta e^{-\delta y}d\delta=1$, we obtain the expression

[TABLE]

The second integral is identical to the third under $\alpha\to 1-\alpha$. It can be computed using an integration by part,

[TABLE]

where we used that $\operatorname{erfc}(x)=\frac{2}{\sqrt{\pi}}\int_{x}^{\infty}dze^{-z^{2}}$. Note that the last term of this equation allows to simplify half of the first term of Eq. (93). As there are two such terms (coming from the integrals with $\alpha$ and $1-\alpha$), the final result reads

[TABLE]

where we used that $\operatorname{erfc}(0)=1$ and $\operatorname{erfc}(x\to\infty)=0$.

Appendix C Computations of $M_{p}$ and $m_{p}$

To obtain the value of $M_{p}$, we compute the moment of order $p$ of the large deviation scaling function of the PDF $\Psi(\tilde{\Delta})$. Using that $\Psi(x)=(\sqrt{2}/\pi)\partial_{x}^{2}(\sinh^{-1}(\sqrt{2}x))$, $M_{p}$ reads after integration by part

[TABLE]

Changing variable from $x\to z=\sqrt{2}x$, we obtain

[TABLE]

Finally, using the integral representation of the Riemann Zeta function [49],

[TABLE]

we obtain the final result in Eq. (79).

To obtain the value of $m_{p}$, we proceed similarly, computing the moment of order $p$ of the large deviation scaling function of the PDF $\varphi_{0}(\tilde{\Delta})$. Using that $\varphi_{0}(x)=(8\pi)^{-1/2}\partial_{x}^{2}(\sinh^{-2}(x))$, $m_{p}$ reads after integration by part

[TABLE]

Changing variable $x\to z=\sqrt{2}x$, we obtain

[TABLE]

Finally, using the integral representation [50] of the Riemann Zeta function,

[TABLE]

we obtain the final expression in Eq. (80).