Extreme value statistics of ergodic Markov processes from first passage times in the large deviation limit

TL;DR

This paper develops a framework linking extreme value statistics of ergodic Markov processes to their relaxation eigenspectra, providing bounds and approximations for long-time extreme value distributions in confining potentials.

Contribution

It introduces a novel approach connecting extreme value functionals to relaxation eigenspectra and derives bounds for first passage densities in ergodic reversible Markov processes.

Findings

Bounds on long-time first passage densities derived

Large deviation limit accurately approximates extreme value statistics

Applicable to Ornstein-Uhlenbeck and Bessel processes

Abstract

Extreme value functionals of stochastic processes are inverse functionals of the first passage time -- a connection that renders their probability distribution functions equivalent. Here, we deepen this link and establish a framework for analyzing extreme value statistics of ergodic reversible Markov processes in confining potentials on the hand of the underlying relaxation eigenspectra. We derive a chain of inequalities, which bounds the long-time asymptotics of first passage densities, and thereby extrema, from above and from below. The bounds involve a time integral of the transition probability density describing the relaxation towards equilibrium. We apply our general results to the analysis of extreme value statistics at long times in the case of Ornstein-Uhlenbeck process and a 3-dimensional Brownian motion confined to a sphere, also known as Bessel process. We find that even on time-scales that are shorter than the equilibration time, the large deviation limit characterizing long-time asymptotics can approximate the statistics of extreme values remarkably well. Our findings provide a novel perspective on the study of extrema beyond the established limit theorems for sequences of independent random variables and for asymmetric diffusion processes beyond a constant drift.

Equations

$${\dot{x}}_t=-\beta D{U}^{\prime }(x_t)+{\xi }_t$$

$$\partial_{t}P(x,t|x_{0})=\partial_{x}\big{[}\partial_{x}+U^{\prime}(x)\big{]}P(x,t|x_{0})\equiv\mathcal{L}_{\rm FP}P(x,t|x_{0}),$$

$${\overline{m}}_t\equiv \underset{0\le \tau \le t}{\max }(x_{\tau })\text{and}{\underline{m}}_t\equiv \underset{0\le \tau \le t}{\min }(x_{\tau }),$$

$$t_a(x_0)\equiv \underset{\tau \ge 0}{\min }(\tau |x_{\tau }=a)$$

$$P^{\max }(a|t,x_0)\equiv \mathrm{Prob}\left[a>{\overline{m}}_t\right]=\mathrm{Prob}\left[t_a(x_0)>t\right]\equiv {\mathcal{P}}_a(t|x_0),$$

$${\mathcal{P}}_a(t|x_0)={\int }_t^{\infty }{\wp }_a(\tau |x_0)\mathrm{d}\tau ,$$

$$p^{\kappa}(a|t,x_{0})=\big{|}\partial_{a}P^{\kappa}(a|t,x_{0})\big{|}=\big{|}\partial_{a}\mathcal{P}_{a}(t|x_{0})\big{|},$$

$$⟨{\overline{m}}_t⟩$$

$$⟨{\underline{m}}_t⟩$$

$$P(a,t|x_0)={\int }_0^{t}P(a,t-\tau |a){\wp }_a(\tau |x_0)\mathrm{d}\tau ,$$

$${\tilde{\wp }}_a(s|x_0)=\frac{\tilde{P}(a,s|x_0)}{\tilde{P}(a,s|a)},$$

$$[{\mathcal{L}}_{\mathrm{FP}}-s]\tilde{P}(x,s|x_0)=-\delta (x-x_0)$$

$$\tilde{P}(a,s|x_0)=\frac{P^{\mathrm{eq}}(a)}{s}+\sum _{k=1}^{\infty }\frac{{\psi }_k(a){\psi }_k(x_0)/P^{\mathrm{eq}}(x_0)}{s+{\lambda }_k},$$

$${\tilde{\wp }}_a(s|x_0)=\sum _{k=1}^{\infty }\frac{w_k(a,x_0){\mu }_k(a)}{{\mu }_k+s}.$$

$$w_{k}(a,x_{0})\mu_{k}(a)=\frac{\tilde{P}(a,s|x_{0})}{\partial_{s}\tilde{P}(a,s|a)}\bigg{|}_{s=-\mu_{k}}.$$

$${\wp }_a(t|x_0)=\sum _{k>0}{w}_k(a,x_0){\mu }_k(a){\text{rme}}^{-{\mu }_k(a)t},$$

$${\mathcal{P}}_a(t|x_0)={\int }_t^{\infty }{\wp }_a(\tau |x_0)\mathrm{d}\tau =\sum _{k>0}{w}_k(a,x_0){\text{rme}}^{-{\mu }_k(a)t},$$

$$p^{\kappa }(a|t,x_0)$$

$$p^{\kappa}_{\rm LD}(a|t,x_{0})\equiv\big{|}\partial_{a}w_{1}(a,x_{0})\rme^{-\mu_{1}(a)t}\big{|}\simeq p^{\kappa}(a|t,x_{0}),$$

$$f(s)=s\tilde{P}(a,s|a)=\sum _{n\ge 0}{\sigma }_n{s}^n/n!$$

$${\mu }_1(a)=\sum _{n=1}^{\infty }\frac{{\sigma }_0^n}{{\sigma }_1^{2n-1}}\frac{\det {\mathcal{A}}_n}{(n-1)}$$

$$\mathcal{A}_{n}^{ij}\equiv\frac{\sigma_{i-j+2}\Theta(i-j+1)}{(i-j+2)!}\times\left\{\begin{array}[]{@{}ll@{}}i&\quad\text{if $j=1$,}\\ n(i-j+1)+j-1&\quad\text{if $j>1$,}\end{array}\right.$$

$$w_{1}(a)=\frac{\tilde{P}(a,-\mu_{1}(a)|x_{0})}{\mu_{1}(a)\partial_{s}\tilde{P}(a,s|a)}\bigg{|}_{s=-\mu_{1}(a)}.$$

$${\tilde{\mu }}_1(a)\equiv \frac{{\sigma }_0}{{\sigma }_1}\approx {\mu }_1(a),$$

$${\sigma }_0=P^{\mathrm{eq}}(a)\text{and}{\sigma }_1={\int }_0^{\infty }[P(a,t|a)-P^{\mathrm{eq}}(a)]\mathrm{d}t.$$

$$\frac{{\tilde{\mu }}_1(a)}{1+{\tilde{\mu }}_1(a)/{\lambda }_1}=\frac{{\lambda }_1{\sigma }_0}{{\sigma }_0+{\lambda }_1{\sigma }_1}\le {\mu }_1(a)\le \frac{{\sigma }_0}{{\sigma }_1}={\tilde{\mu }}_1(a),$$

$$P(a,t|x_{0})=\frac{1}{\sqrt{2\pi(1-\rme^{-2t})}}\exp\bigg{[}-\frac{(a-x_{0}\rme^{-t})^{2}}{2(1-\rme^{-2t})}\bigg{]}$$

$$\tilde{\mu}_{1}(a)=\int_{0}^{\infty}\bigg{[}\frac{1}{\sqrt{1-\rme^{-2t}}}\exp\bigg{(}\frac{a^{2}\rme^{-t}}{1+\rme^{-t}}\bigg{)}-1\bigg{]}{\rm d}t\simeq\mu_{1}(a).$$

$${\tilde{\mu }}_1(1+{\tilde{\mu }}_1{)}^{-1}\le {\mu }_1\le {\tilde{\mu }}_1,$$

$$\tilde{P}(a,s|x_0)=\Gamma (s)2^s{P}^{\mathrm{eq}}(a)H_{-s}(-a/\sqrt{2})H_{-s}(x_0/\sqrt{2}),$$

$${\tilde{\wp }}_a(s|x_0)=\frac{H_{-s}(-x_0/\sqrt{2})}{H_{-s}(-a/\sqrt{2})},$$

Full text

david.hartich@mpibpc.mpg.de,agodec@mpibpc.mpg.de

1 Introduction

The statistical properties of extreme values, which correspond to record-breaking events of a stochastic process, attracted increasing interest in various fields of research over the past decades. For example, climate changes were found to be reflected in the appearance of extreme (record-breaking) temperatures [1, 2], rainfall [3, 4], and possibly other extreme weather conditions [5]. Statistics of records are also important in the context of earthquakes [6], in studies of stock pricing in economics [7, 8], sports [9, 10], and in the theory of random matrices [11, 12, 13] to name but a few (see, e.g., Ref. [14, 15] for a more detailed overview).

In sequences of independent random variables extreme values approach one of the three classes of limiting distributions, which are denoted by the Gumbel [16], Fréchet, and Weibull distributions (see, e.g., Refs. [17, 18, 19, 20, 14]). However, as soon as consecutive time steps of a stochastic process become correlated, a theoretical discussion of the statistics of extrema becomes more challenging [21]. In this case universal laws have been discovered, for example, for processes with symmetric step-length distributions [22, 23]. Subsequent studies also investigated correlations between records [24, 25] as well as their persistence [26, 27, 28] and number [29], and extensions have been made to processes with constant drift [30, 27, 31, 32] (see also Ref. [33] for an interesting experiment with trapped Cs atoms). A recent physical application includes the observation that the mean value of the minimum of the entropy production in stationary driven systems is bounded by the negative value of Boltzmann’s constant “$-k_{\rm B}$” [34], which is also confirmed by experiments [35].

More broadly, a deep and important connection has been established, relating the statistics of extreme values to first passage times [23, 36, 14, 37]. In this work we deepen this connection between the first passage and the extremum functional, which allows us to obtain the statistics of extreme values in finite time for Markovian diffusion processes in confining potentials on time-scales, where consecutive time-steps remain correlated. Exploiting further a duality between first passage processes and ensemble propagation [38, 39] we derive a chain of inequalities, which bound the long time asymptotics (i.e., the large deviation limit) of the probability densities of extrema both from above and from below. As we will show, the large deviation limit approximates the probability density of extreme values surprisingly well even on relatively short time-scales.

The paper is organized as follows. In Sec. 2 we recapitulate the well-known connection between distributions of extrema and first passage time densities. We then utilize recent findings on the large time asymptotics of first passage time densities [38, 39] to determine distributions of extrema in the large deviation limit. The usefulness of our general results is demonstrated in Sec. 3, by determining the long-time statistics of maxima of the Ornstein-Uhlenbeck process and the statistics of the minimum of the 3-dimensional Brownian motion (Bessel process) confined to a sphere. All analytical results are corroborated by Brownian dynamics simulations. We conclude in Sec. 4.

2 Fundamentals

2.1 Extreme values from first passage times

We consider processes governed by an overdamped Langevin equation

[TABLE]

where $U^{\prime}(x)=\partial_{x}U(x)$ is the gradient of a potential $U(x)$ and $\xi_{t}$ stands for Gaussian white noise with zero mean and covariance $\langle\xi_{t}\xi_{t^{\prime}}\rangle=2D\delta(t-t^{\prime})$. Without any loss of generality we set the inverse temperature $\beta$ and diffusion coefficient $D$ to unity ($\beta\equiv D\equiv 1$), i.e., free energies $U$ are expressed in units of $k_{\rm B}T$. The Fokker-Planck equation corresponding to the Langevin equation (1) reads [40]

[TABLE]

where $P(x,t|x_{0})=\langle\delta(x-x_{t})\rangle$ is the normalized probability density for a particle starting from $x_{0}$ to be found at position $x_{t}=x$ at time $t$ with the initial condition $P(x,0|x_{0})=\delta(x-x_{0})$. The probability density function relaxes to the normalized Boltzmann-Gibbs equilibrium density $P^{\rm eq}(x)\equiv P(x,\infty|x_{0})\propto\rme^{-U(x)}$ for any $x_{0}$, which requires a sufficiently confining potential $U(x)$.

While the probability density $P(x,t|x_{0})$ only depends on the initial and final states, the extreme values are functionals that depend on the entire history along a trajectory $\{x_{\tau}\}_{0\leq\tau\leq t}$. We define the maximum and the minimum of the process $x_{t}$ as

[TABLE]

respectively. For a given initial condition $x_{0}$ the extrema satisfy $\overline{m}_{t}\geq x_{0}$ as well as $\underline{m}_{t}\leq x_{0}$, where $\overline{m}_{t}$ is non-decreasing and $\underline{m}_{t}$ non-increasing in time $t$ (see Fig. 1). It can be shown that the first passage time defined as

[TABLE]

is an inverse functional of the extrema. To see that we consider the distribution function of the maximum of the process $\overline{m}_{t}=\max_{0\leq\tau\leq t}(x_{\tau})$, i.e., the probability that $\overline{m}_{t}$ exceeds the value $a\geq x_{0}$, $P^{\rm max}(a|t,x_{0})$, which satisfies

[TABLE]

where $\mathcal{P}_{a}(t|x_{0})$ is called the survival probability (see also Fig. 1a).

Eq. (5) can be interpreted as follows: each path, whose maximum $\overline{m}_{t}$ after time $t$ is smaller than $a$, must have a first passage time from $x_{0}$ to $a\geq x_{0}$, $t_{a}(x_{0})$, larger than $t$. Eq. (5) connects the first passage time functional $t_{a}(x_{0})$ (where time is stochastic and position is fixed to $a$) to the maximum value functional $\overline{m}_{t}$ (where time is fixed to $t$ and the position is stochastic). Note that the survival probability can also be expressed as the integral over the first passage time density $\wp_{a}$ via

[TABLE]

i.e., $\wp_{a}(\tau|x_{0})=-\partial_{t}\mathcal{P}_{a}(t|x_{0})$.

The minimum of a process $\underline{m}_{t}=\min_{0\leq\tau\leq t}x_{t}$ can be studied in a similar manner as the maximum, since $\underline{m}_{t}$ is equivalent to the maximum of the reflected process $-x_{\tau}$, i.e., $\underline{m}_{t}=-\max_{0\leq\tau\leq t}(-x_{\tau})$. Hence, in the case of the minimum ($a\leq x_{0}$) Eq. (5) holds with the replacement $P^{\rm min}(a|t,x_{0})={\rm Prob}\left[a<\underline{m}_{t}\right]=\mathcal{P}_{a}(t|x_{0})$, which is illustrated in Fig. 1b. For convenience, we simply refer to $P^{\kappa}$ as the extremum distribution function, which in the case $x_{0}\leq a$ corresponds to the maximum distribution ($\kappa={\rm max}$) and for $x_{0}\geq a$ to the minimum distribution ($\kappa={\rm min}$). The density of the extremum ($\kappa={\rm max},{\rm min}$) in either case is then given by the slope of the distribution function

[TABLE]

where in the second step we used Eq. (5). Eq. (7) describes the density of maxima ($\kappa={\rm max}$) for $a\geq x_{0}$ and the density of the minima ($\kappa={\rm min}$) for $a\leq x_{0}$. For example, the mean value of the maximum and minimum are given, respectively, by

[TABLE]

where we used Eq. (7) and performed a partial integration in the last step in both lines.

In the following subsection we focus on the probability density function of the two extrema $p^{\rm max}$ and $p^{\rm min}$, whereas further discussions on the mean of extreme value fluctuations ($\langle\underline{m}_{t}\rangle$ or $\langle\overline{m}_{t}\rangle$) can be found, for example, in Refs. [41, 21, 23, 24, 25]. Notably, $t_{a}(x_{0})$ – the first passage time from $x_{0}$ to $a\geq x_{0}$ – is unaffected by the potential landscape $U(x)$ beyond $x\geq a$, which according to Eq. (5) implies that any two potentials $U_{1}(x),U_{2}(x)$ with $U_{1}(x)=U_{2}(x)$ for all $x\leq R$ generate the same maximum distribution $P^{\rm max}$ for all $a\leq R$.

2.2 First passage time statistics from ensemble propagation

According to Eq. (5) the problem of determining the statistics of the extremum $P^{\rm max}(a|t,x_{0})$ is in fact equivalent to determining the survival probability $\mathcal{P}_{a}(t|x_{0})$, or, according to Eq. (6), to determining the first passage time density $\wp_{a}(t|x_{0})=-\partial_{t}\mathcal{P}_{a}(t|x_{0})$, which will be the central goal of this section.

We determine the first passage time density (or survival probability) using the renewal theorem [42]

[TABLE]

reflecting the fact that all the paths starting from $x_{0}$ and ending up in $a$ after time $t$ by construction must reach $a$ for the first time at some time $\tau\leq t$, and then return to $a$ again after time $t-\tau$. We have recently established a duality between first passage and relaxation processes, i.e., between $P$ and $\wp_{a}$, which will allow us to solve Eq. (9) for $\wp_{a}$ in the following manner (see also Refs. [38, 39] for more details).

First, we Laplace-transform111The Laplace transform of a function $f(t)$ is defined by $\tilde{f}(s)=\int_{0}^{\infty}f(t)\rme^{-st}{\rm d}t$. the renewal theorem (9) in time $(t\to s)$, which converts the convolution to a product, $\tilde{P}(a,s|x_{0})=\tilde{P}(a,s|a)\tilde{\wp}_{a}(s|x_{0})$, implying

[TABLE]

where $\tilde{\wp}_{a}(s|x_{0})$ is the Laplace transform of the first passage time density and $\tilde{P}(x,s|x_{0})$ obeys the Laplace transformed Fokker-Planck equation (2)

[TABLE]

with natural boundary conditions. The next step is to render Eq. (10) explicit in the time domain, i.e., to find the explicit inverse Laplace transform $\tilde{\wp}_{a}(s|x_{0})\to\wp_{a}(t|x_{0})$.

Therefore, recalling that we consider sufficiently confining potentials $U(x)$, there exist a spectral expansion of the Fokker-Planck operator $\mathcal{L}_{\rm FP}$ with discrete eigenvalues $-\lambda_{k}\leq 0$ for $k=0,1,\ldots$ and corresponding symmetrized nontrivial solutions $\psi_{k}(x)$ to the eigenequation $\mathcal{L}_{\rm FP}\psi_{k}(x)=-\lambda_{k}\psi_{k}(x)$, which are assumed to be normalized $\int\psi_{k}(x)^{2}/P^{\rm eq}(x){\rm d}x\equiv 1$. The ground state $\psi_{0}$ corresponding to eigenvalue $\lambda_{0}=0$ represents the equilibrium Boltzmann distribution $P^{\rm eq}(x)=\psi_{0}(x)$. Using the eigenfunctions $\{\psi_{k}\}$ and eigenvalues $\{\lambda_{k}\}$ defined this way, the Laplace transform of the ensemble propagator can be written in the form

[TABLE]

where $\psi_{k}(x)/P^{\rm eq}(x)\equiv\psi_{k}^{\dagger}(x)$ are in fact the eigenfunctions to the adjoint of $\mathcal{L}_{\rm FP}$, that is $\mathcal{L}_{\rm FP}^{\dagger}\psi^{\dagger}(x)=-\lambda_{k}\psi^{\dagger}(x)$. Note that in our previous work we used the equivalent non-symmetric eigenspectrum with right and left eigenfunctions $\psi_{k}$ and $\psi_{k}^{\dagger}$, respectively [38, 39, 43]. Since the Laplace transform of a function $f$ with a simple pole $\tilde{f}(s)=(s+\lambda)^{-1}$ yields in the time domain an exponentially decaying function $f(t)=\rme^{-\lambda t}$ with rate $\lambda$, we can interpret the eigenvalues $\lambda_{k}$ as relaxation rates, which characterize the speed at which the dynamics governed by Eq. (2) approaches the equilibrium $P^{\rm eq}(x)\propto\rme^{-U(x)}$.

The Laplace transform of the first passage time density $\tilde{\wp}_{a}(s|x_{0})$, as well has simple poles, which are located at $s=-\mu_{k}$ ($k=1,2,\ldots$) and need to be determined for Eq. (10) to be written as [38]

[TABLE]

The expansion in Eq. (13) can formally be found by determining the zeros $s=-\mu_{k}$ that solve $\tilde{P}(a,s|a)=0$, to which we refer as first passage rates $\mu_{k}$. Determining all first passage rates $\mu_{k}$, while doable in general, is rather involved and is described in [38, 39], whereas detailed information on the determination of slowest rate $\mu_{1}$, to which we refer to as large deviation limit, can be found in Sec. 2.3 below as well as in [44]. If all first passage rates $\{\mu_{k}\}$ are known, we can obtain the corresponding weights $w_{k}(a,x_{0})$ in Eq. (13) directly from Eq. (10) using Cauchy’s residue theorem

[TABLE]

Dividing Eq. (14) by $\mu_{k}(a)$ yields the “weights” $w_{k}(a,x_{0})$, which according to Eq. (13) are normalized such that $\tilde{\wp}_{a}(0|x_{0})=\sum_{k}w_{k}(a,x_{0})=1$. Eq. (13) in turn immediately yields the first passage time density

[TABLE]

and the corresponding survival probability

[TABLE]

where we used Eqs. (6) and (15). Inserting Eq. (16) into (7) allows us to rewrite the probability density of the extremum to have value $a$ at time $t$ as

[TABLE]

where $\kappa={\rm max}$ or $\kappa={\rm min}$.

2.3 Large deviation limit

At long times the extremum ($\overline{m}_{t}$ or $\underline{m}_{t}$) will be dominated by extreme fluctuations of the process $x_{t}$ that are not reflected by the “typical” equilibrium measure given by $P^{\rm eq}(x)\propto\rme^{-U(x)}$. As a result, the extreme value distribution may differ substantially from the equilibrium Boltzmann distribution $P^{\rm eq}(x)$. Fortunately, at long times the first passage distribution will be dominated solely by the slowest first passage time-scale $1/\mu_{1}(a)$, which leads to what we refer here to as the large deviation limit that reads

[TABLE]

where “$\simeq$” denotes the asymptotic equality in the limit $t\to\infty$ and $\kappa={\rm max},{\rm min}$. An explicit general method to determine $w_{1}(a,x_{0})$ and $\mu_{1}(a)$ can be found in [44, 38, 39]. We note that the large deviation limit becomes exact in the long time limit $\rme^{-\mu_{1}(a)t}\gg\rme^{-\mu_{2}(a)t}$ as well as whenever $w_{1}(a,x_{0})\gg|w_{k\geq 2}(a,x_{0})|$ holds.

We recall that according to Eq. (13) each first passage rate, $\mu_{k}$, is located at a simple pole (at $s=-\mu_{k}$) of $\tilde{\wp}_{a}(s|x_{0})$, which according to Eq. (10), is also a root of $\tilde{P}(a,s|a)$. Hence, the large deviation limit “$\mu_{1}$” is characterized by the root ($s<0$) closest to the origin, $s=-\mu_{1}$, solving $\tilde{P}(a,s|a)=0$. In order to determine $\mu_{1}$ exactly, we Taylor-expand the function

[TABLE]

around $s=0$, where $\sigma_{n}$ is the $n$th derivative of $f$ with respect to $s$, which according to Eq. (12) holds for all $s$ within the radius of convergence $|s|<\lambda_{1}$.222The radius of convergence is limited by the pole of $f(s)=s\tilde{P}(a,s|a)$ which is closest to the origin. According to Eq. (12) the closest pole to $s=0$ is located at $s=-\lambda_{1}$, yielding a converging sum Eq. (19) for all $|s|<\lambda_{1}$. Note that $f(s)=s\tilde{P}(a,s|a)$ has the same non-trivial roots but, in contrast to $\tilde{P}(a,s|a)$, does not have a pole at the origin (see Eq. (12)), which is why we are always allowed to expand $f$ as in Eq. (19). The closest non-trivial zero $s=-\mu_{1}$ with $0<\mu_{1}\leq\lambda_{1}$ can then be formally be found by a Newton’s iteration, which in terms of a series of almost triangular matrices reads explicitly [39] (see also [44, 38])

[TABLE]

where $\boldsymbol{\mathcal{A}}_{n}$ is a $(n-1)\times(n-1)$ almost triangular matrix with elements ($i,j=1,\ldots,n-1$)

[TABLE]

with $\Theta(l)=1$ if $l\geq 0$ and $\Theta(l)=0$ if $l<0$ as well as $\det\boldsymbol{\mathcal{A}}_{1}\equiv 1$. Eq. (20) exactly determines the first non-trivial root, $f(s)=0$ with $s=-\mu_{1}$, at which the right hand side of Eq. (19) vanishes. It should be noted that determining $\mu_{1}$ in Eq. (20) requires only $\tilde{P}(a,s|a)$ or the coefficients $\sigma_{n}$ from Eq. (19), whereas the expansion Eq. (12) including the eigenvalues $\{\lambda_{k}\}$ is generally not required to be known. The weight $w_{1}(a)$ can then be deduced from Cauchy’s residue theorem Eq. (14)

[TABLE]

Equation (18) with Eqs. (20) and (22) fully characterize the large deviation limit of the density of the extreme value $p^{\rm max}_{\rm LD}(a|x_{0})$.

2.4 Large deviation limit in the presence of a spectral gap

In the large time limit the probability mass of the extremum $p^{\rm max}(a|t,x_{0})$ or $p^{\rm min}(a|t,x_{0})$ concentrates at the potential boundaries (i.e., $U(a)\gg k_{\rm B}T$), such that we can accurately approximate $\mu_{1}$ by truncating Eq. (20) already after the first term yielding (see Ref. [38] for more details)

[TABLE]

where using Eq. (12) we can identify

[TABLE]

Since $f(s)\equiv s\tilde{P}(a,s|x_{0})=\sigma_{0}+\sigma_{1}s+\mathrm{O}(s)^{2}$, we expect Eq. (23) to be quite accurate as soon as the formal condition $\tilde{\mu}_{1}\ll\lambda_{1}$ is met, where $\lambda_{1}$ from Eq. (12) is the slowest rate at which the system approaches the equilibrium density [38]. Note that $\lambda_{1}$ in fact does not need to be known, as Eq. (23) necessarily becomes accurate at sufficiently high potential values $U(a)$, such that $a$ is not located in the deepest point in the potential [38]. This also follows from the work of Matkowsky and Schuss, who have shown that $\lambda_{1}$ is the expected time to overcome the barriers on the way to the deepest potential well [45].

In fact, at very long times $t\to\infty$ the probability mass $p^{\kappa}(a|t,x_{0})\simeq p^{\kappa}_{\rm LD}(a|t,x_{0})$ (with $\kappa={\rm\max},{\rm min}$) will inevitably be pushed towards the boundaries with high potential values, which will again render Eq. (23) asymptotically exact in the limit $U(a)\to\infty$. To prove that Eq. (23) indeed becomes asymptotically exact, we inspect Eq. (12) in the following way. First, we find that $f(s)=s\tilde{P}(a,s|a)$ is a concave function $f^{\prime\prime}(s)\leq 0$ within the interval $-\lambda_{1}\leq s\leq 0$, whereas $g(s)\equiv s(s+\lambda_{1})\tilde{P}(a,s|a)=(s+\lambda_{1})f(s)$ is a convex function $g^{\prime\prime}(s)\geq 0$ in the same interval. Hence, we find that the tangent $t_{f}(s)=\sigma_{0}+s\sigma_{1}$ to $f$ and the tangent $t_{g}(s)=\lambda_{1}\sigma_{0}+(\sigma_{0}+\lambda_{1}\sigma_{1})s$ to $g$ have roots that sandwich $s=-\mu_{1}$ according to

[TABLE]

where the lower bound, $s=-\tilde{\mu}_{1}(1+\tilde{\mu}_{1}/\lambda_{1})^{-1}$, solves $t_{g}(s)=0$ and the upper bound, $s=-\tilde{\mu}_{1}$, solves $t_{f}(s)=0$. The chain of inequalities Eq. (25) and its implications, which we explore below, are the main result of this paper. Notably, in the limit of $U(a)\gg k_{\rm B}T$, where $\tilde{\mu}_{1}\to 0$ (i.e., $\tilde{\mu}_{1}\ll\lambda_{1}$ [45]) holds, the inequalities in Eq. (25) saturate and provide an asymptotically exact value for $\mu_{1}$.

We emphasize that the chain of inequalities (25) holds for the slowest time-scale $\mu_{1}^{-1}$ of the first passage process as well as for the slowest time-scale of the extremum functional, i.e. either the maximum or the minimum. For example, if $\mu_{1}^{\text{min}}(a)$ and $\mu_{1}^{\text{max}}(a)$ denote the large deviation limit of the minimum and maximum, respectively, then the slowest first passage rate is given by $\mu_{1}(a)=\min[\mu_{1}^{\text{max}}(a),\mu_{1}^{\text{min}}(a)]$. Since the maximum $\overline{m}_{t}=a\geq x_{0}$ after a long time $t$ will be more likely located at the “right” border of a confining potential, where $U(\overline{m}_{t})\gg U(x_{0})$, whereas the minimum $\underline{m}_{t}=a\leq x_{0}$ will more likely move to the “left” border, where $U(\underline{m}_{t})\gg U(x_{0})$, we will use Eq. (25) to determine the minimum near the left boundary $a<x_{0}$ and to determine the maximum if $a>x_{0}$ is closer to the right boundary.

To be more specific, we use $\tilde{\mu}_{1}(a)\simeq\mu_{1}^{\text{max}}(a)$ self-consistently for the large deviation limit of the maximum, whenever $\tilde{\mu}_{1}(a)\searrow$ is monotonically decreasing with increasing $a\nearrow$, whereas we use $\tilde{\mu}_{1}(a)\simeq\mu_{1}^{\text{min}}(a)$ for the large deviation limit of the minimum, whenever $\tilde{\mu}_{1}(a)\searrow$ is monotonically decreasing with decreasing $a\searrow$. Notably, if $a$ is located at a reflecting boundary, where formally $U(x)=\infty$ for $x\leq a$, we immediately get $\mu_{1}(a)=\mu_{1}^{\text{min}}(a)$ and $\mu_{1}^{\text{max}}(a)=\infty$, since for any $x_{0}>a$ the maximum $\overline{m}_{t}$ cannot reach any value below $x_{0}$ and hence $a$ certainly cannot correspond to the maximum.333A similar finding can be found in [39], where $\mu_{1}(a)=\mu_{1}^{\text{min}}(a)$ corresponds to a first passage to $a$, entering from the right, and $\mu_{1}(a)=\mu_{1}^{\text{max}}(a)$ corresponds to a first passage to $a$, entering from the left. For example, if a reflecting boundary is located at $a$ with $U(x)=\infty$ for $x\leq a$, it is impossible to enter $a$ from the left.

3 Examples

3.1 Statistics of maxima in the Ornstein Uhlenbeck process

As our first example we consider the Ornstein-Uhlenbeck process with $U(x)=x^{2}/2$. The corresponding propagator in the time domain is well known and reads [40]

[TABLE]

with a Gaussian equilibrium density $P^{\rm eq}(a)=P(a,\infty|x_{0})=(2\pi)^{-1/2}\exp(-a^{2}/2)$. Inserting Eq. (26) into Eq. (23) yields for $a\geq 0$ the following approximation for the large deviation eigenvalue of the density of the maximum444For $a\leq 0$ Eq. (27) approximates the large deviation limit of the minimum functional.

[TABLE]

The relaxation eigenvalues are integers $\lambda_{k}=k$ with $k=0,1,\dots$, such that Eq. (25) translates into

[TABLE]

where the upper limit $\tilde{\mu}_{1}$ is depicted in Fig. 2a as the dash-dotted red line, the lower limit $\tilde{\mu}_{1}(1+\tilde{\mu}_{1})^{-1}$ as the dashed green line, and the exact value $\mu_{1}$, determined as described below, is given by the solid blue line. The inset displays the same results but scaled by the exact value $\mu_{1}$. We emphasize that it is not necessary to determine $\mu_{1}$ in order to show that Eq. (27) asymptotically saturates when $\tilde{\mu}_{1}$ approaches zero in the limit of large $a$, since $\tilde{\mu}_{1}\simeq\mu_{1}$ follows immediately from $\tilde{\mu}_{1}(1+\tilde{\mu}_{1})^{-1}\simeq\tilde{\mu}_{1}$ (for $\tilde{\mu}_{1}\to 0$) as well as from Eq. (28). For completeness, we also present in Fig. 2a (dotted line) the long time asymptotics, $\mu_{1}\simeq(2\pi)^{-1/2}a\rme^{-a^{2}/2}$, for the limit $a\to\infty$, which have been reported previously [46, 47, 44].

In order to determine the large deviation eigenvalue $\mu_{1}$ and weight $w_{1}$ entering $p^{\rm max}_{\rm LD}(a|t,x_{0})=\partial_{a}w_{1}(a,x_{0})\rme^{-\mu_{1}(a)t}$, we Laplace-transform the propagator (26) in time ($t\to s$), which for $x_{0}\leq a$ yields (see also [42])

[TABLE]

where $\Gamma(s)$ is the complex gamma function and $H_{s}(y)$ is the generalized Hermite polynomial. Inserting Eq. (29) into the renewal theorem (10) yields [42, 49, 50]

[TABLE]

where $s=-\mu_{1}$ is the root, $H_{-s}(-a/\sqrt{2})=0$, closest to the origin such that the weight in Eq. (22) becomes [46, 48]

[TABLE]

where we introduced $h_{s}(y)\equiv\partial_{s}H_{s}(y)$. For convenience we determined $\mu_{1}$ and $w_{1}$ numerically according to Ref. [48].555We note that with the eigenfunctions $\psi_{k}(a)=P^{\rm eq}(a)(k!2^{k})^{-1/2}H_{k}(a/\sqrt{2})$ (see, e.g., Ref. [40]) and Eqs. (12) and (19) we can formally identify

$\sigma_{0}=P^{\rm eq}(a)\qquad\text{and}\qquad\frac{\sigma_{n}}{n!}=P^{\rm eq}(a)(-1)^{n+1}\sum_{k=1}^{\infty}\frac{H_{k}(a/\sqrt{2})^{2}}{k!2^{k}k^{n}},$

which with Eq. (20) would be an alternative but equivalent approach for determining $\mu_{1}$ as done, e.g. in Ref. [39].

The results in Fig. 2a confirm the validity of the chain of inequalities in Eq. (25), which, as already mentioned, become asymptotically tight in the limit of high values of the potential, $U(a)=a^{2}/2\gg 1$.

The lines in Fig. 2b represent the large deviation limit of the density of the maximum $p^{\rm max}_{\rm LD}(a|x_{0})=\partial_{a}w_{1}(a,x_{0})\rme^{-\mu_{1}(a)t}$, which agree rather well with the density of the maximum $p^{\rm max}(a|x_{0})$ (symbols) obtained from Brownian dynamic simulations with a time step ${\rm d}t=10^{-5}$ using $10^{5}$ trajectories.

We note that the large deviation limit $p^{\rm max}_{\rm LD}$ (see lines in Fig. 2b) approximates the density of the maximum $p^{\rm max}$ quite well already on relatively short time-scales $t\sim\lambda_{1}^{-1}=1$ (see triangles and dash-dotted light green line), where $\lambda_{1}^{-1}$ represents the equilibration time of the Ornstein-Uhlenbeck process. Notably, even for long times (see, e.g., $t=50$ in Fig. 2b), the left and right tails of the density of the maximum remain asymmetric yet still deviating from a Gumbel distribution [16] (see thick gray line in Fig. 2b). This indicates that the extreme value theorem for sequences of uncorrelated random variables becomes valid on much longer time-scales. Therefore, the large deviation limit presented here allows us to approximate extreme value statistics exceptionally well despite the fact that the extreme value theorem does not yet apply.666 We find that the probability density of the maximum approaches a Gumbel density on extremely large time-scales $t\gtrsim 10^{3}$. The underlying assumptions are the approximation $\mu_{1}\simeq a(2\pi)^{-1/2}\rme^{-a^{2}/2}$ (see inset of Fig. 2a for deviations) and $w_{1}\simeq 1$ (which holds for $a\gtrsim 3$).

3.2 Density of the minimum of the confined Bessel process

In our second example we consider the minimal distance to the origin of Brownian motion inside a $d$-dimensional sphere with inner radius $R_{-}\geq 0$ and a reflecting boundary at $R_{+}<\infty$ (see Fig. 3a for an illustration with $R_{-}=0$ and $d=2$). The distance from the origin $x_{t}$ (i.e. the radius) at time $t$ within the interval $R_{-}\leq x\leq R_{+}$ obeys the Langevin equation

[TABLE]

where $\langle\xi_{t}\rangle$ and $\langle\xi_{t}\xi_{t^{\prime}}\rangle=2\delta(t-t^{\prime})$. This process is also known as the Bessel process [51, 52]. We note that the maximum excursion of the free Bessel process (see e.g. [53]), which in the present context corresponds to the limiting case with $R_{+}=\infty$ and will not be considered here, allows in the specific case of $d=3$ a mapping onto a simpler problem for the 1-dimensional Brownian motion [51].

Comparing Eq. (1) and Eq. (32) allows us to identify the geometric free energy $U(x)=-(d-1)\ln x$ of purely entropic origin and accounts for the invariance with respect to angular degrees of freedom $\propto x^{d-1}$ (see Fig. 3b). The equilibrium measure corresponds to a uniform distribution in a $d$-dimensional hyperspherical shell and is given by $P^{\rm eq}(x)=dx^{d-1}/(R_{+}^{d}-R_{-}^{d})$.

For simplicity we here from restrict our discussion to the case $d=3$, yielding the Fokker-Planck equation

[TABLE]

with zero flux boundary condition $J(R_{\pm},t|x_{0})=0$, where $J(x,t|x_{0})\equiv(2/x-\partial_{x})P(x,t|x_{0})$. We emphasize that the probability density $P$ is normalized according to $\int_{R_{-}}^{R_{+}}P(x,t|x_{0}){\rm d}x=1$, whereas the radial density, discussed for example in [54], would correspond to $P(x,t|x_{0})/(4\pi x^{2})$ instead. A Laplace transform in $t$ yields

[TABLE]

where the solution $\tilde{P}(x,s|x_{0})$ can be constructed from the two solutions of the homogeneous problem, $v_{1}(x,s)=x\rme^{-x\sqrt{s}}$ and $v_{2}(x,s)=x\rme^{x\sqrt{s}}/\sqrt{s}$ obtained by setting the right hand side of Eq. (34) to zero. The Laplace transform of the propagator for a Brownian particle confined between $R_{-}=a$ and $R_{+}=R$ in turn reads

[TABLE]

We are allowed to choose $R_{-}=a$, since the first passage time distribution from $x_{0}$ to $a\leq x_{0}$ is not affected by the potential $U(x)$ in the region $x\leq a$, where $R_{-}=a$ formally corresponds to $U(x)=\infty$ for $x\leq a$ (see also the discussion at the end of Sec. 2.1). Most importantly, setting $R_{-}=a$ removes all roots of $\tilde{P}(a,s|a)$, which would account for the maximum of the Bessel process. In other words, in the presence of a reflecting boundary at $a$ every single root of $\tilde{P}(a,s|x_{0}{=}a)$ from Eq. (35) is indeed a first passage time scale for approaching $a$ for the first time from above (for more details on the influence of boundary condition please see [39]). Moreover, the limit $R=R_{+}=\infty$, which is not considered here, would allow us to map the 3d-Bessel process to 1d Brownian motion [51] with $\tilde{\wp}_{a}(s|x_{0})=\tilde{P}(a,s|x_{0})/\tilde{P}(a,s|a)\to(a/x_{0})\exp[\sqrt{s}(a-x_{0})]$, which would in turn yield the Levy-Smirnov density777The Levy-Smirnov density is defined as $\wp_{a}(t|x_{0})=(a/x_{0})\times(x_{0}-a)/\sqrt{4\pi t^{3}}\times\rme^{-(x_{0}-a)^{2}/(4t)}$..

For $R<\infty$, we use Eq. (35) to identify the Taylor coefficients of $s\tilde{P}(a,s|a)$, which we denote by $\sigma_{n}$ according to Eq. (19). The exact smallest eigenvalue $\mu_{1}$ is then determined using Eq. (20). The results are presented in Fig. 4a (see solid blue line), where we also compare $\mu_{1}$ to the approximation $\tilde{\mu}_{1}$ from Eq. (23) (see dash-dotted red line), which for the 3$d$-Bessel process reads

[TABLE]

Eq. (36) delivers the exact value for $\mu_{1}$ in the limit $a\to 0$ as shown in Fig. 4a, where the relative deviation $\epsilon\equiv(\tilde{\mu}_{1}-\mu_{1})/\mu_{1}$ vanishes in the limit $a\to 0$ (see dashed green line). It should be noted that $\mu_{1}$, given by the series Eq. (20), is in fact an explicit solution of the transcendental equation $R\sqrt{\mu_{1}}=\tan[(R-a)\sqrt{\mu_{1}}]$ in the form of a Newton’s series.

To rationalize why Eq. (36) becomes asymptotically exact as $a$ approaches zero, we recall that $\lambda_{1}\approx 0.81$ is the slowest relaxation rate corresponding to $R_{-}=0$, which solves $R\sqrt{\lambda_{1}}=\tan(R\sqrt{\lambda_{1}})$ (here using $R=5$). Since Eq. (35) obeys a reflecting boundary condition at $a$, we have that $\lambda_{1}(a)\geq\lambda_{1}\approx 0.81$, i.e. the eigenvalue $\mu_{1}(a)$ is bounded by $\tilde{\mu}_{1}(1-\tilde{\mu}_{1}/0.81)\leq\mu_{1}\leq\tilde{\mu}_{1}$. Hence for asymptotically high potentials (here $U(a)=-2\ln a\to\infty$ as $a\to 0$) the inequality Eq. (25) renders $\tilde{\mu}_{1}$ asymptotically exact as soon as $\tilde{\mu}_{1}(a)/0.81\to 0$.

In Fig. 4b we compare the density of the minimum $p^{\rm min}$ (symbols), obtained from simulations of $10^{5}$ trajectories (with $R_{-}=0$, $R_{+}=R=5$ and starting condition $x_{0}=3$), to the corresponding large deviation limit $p^{\rm min}_{\rm LD}(a|t,x_{0})=-\partial_{a}w_{1}(a,x_{0})\rme^{-\mu_{1}(a)t}$ (lines), where we determined $w_{1}$ using Eq. (22) and took the exact $\mu_{1}(a)$ obtained using Eq. (35). The error bars in the simulation results denote 95 % confidence intervals.

By design the large deviation limit $p^{\rm min}_{\rm LD}$ approaches the density of the minimum $p^{\rm min}$ in the long time limit, which is perfectly corroborated by simulation results for $t=40\,(\gg\lambda_{1}^{-1}\approx 1.2)$ in Fig. 4b. Notably, $p^{\rm min}_{\rm LD}$ (see solid dark-blue line in Fig. 4b) approximates quite well the full probability density of the minimum $p^{\rm min}$ (see filled circles).

To our surprise, the large deviation limit $p^{\rm min}_{\rm LD}$ can approximate $p^{\rm min}$ even for smaller times (e.g., $t=1$), i.e. those that are shorter than the equilibration time $\lambda_{1}^{-1}\approx 1.2$ (cf. dash-dotted green line vs. open triangles), which can be explained as follows. For any $a$ within $0\leq a\leq 3=x_{0}$ we find a spectral gap $\mu_{1}(a)\ll\mu_{2}(a)$, which for $1.5\lesssim a\leq 3$ also satisfies $\mu_{2}(a)\gg\lambda_{1}(a)\approx 0.81$. This in turn implies that for $t=1\sim\lambda_{1}^{-1}$ the condition $\rme^{-\mu_{1}t}\gg\rme^{-\mu_{2}t}$ is still met, whereas the relative deviations between $p^{\rm min}$ and $p^{\rm min}_{\rm LD}$ (i.e., between the open triangles and the dash-dotted green line) become substantial for small values of $a$ (see $a\leq 1$) and $\mu_{2}(a)$ concurrently approaches $\lambda_{1}\simeq 0.8$. Once the time exceeds $\lambda_{1}t\approx 0.81t\gg 1$, the condition $\rme^{-\mu_{1}t}\gg\rme^{-\mu_{2}t}$ is satisfied for any value of $a$, and thus $p^{\rm min}_{\rm LD}$ approximates $p^{\rm min}$ over the full range (see symbols and lines in Fig. 4b for $t\geq 5$). Therefore, $p^{\rm min}_{\rm LD}$ approximates $p^{\rm min}$ rather well for any value $a$ and on all time scales longer than the equilibration time scale ($t\gg\lambda_{1}^{-1}$).

Let us finally discuss the “ultimate” long time limit ($t\to\infty$) in which the density of the minimum $p^{\rm min}(a|t,x_{0})$ will be sharply peaked around the shortest distance $a=0$. Inspecting Eq. (36) one can easily find $\mu_{1}(a)=3a/R^{3}+\mathrm{O}(a)$. Moreover, at high values of the potential the weight becomes $w_{1}(a,x_{0})=1+\mathrm{O}(a)$ implying that the limiting density becomes $p^{\rm min}(a|t,x_{0})\to 3tR^{-3}\rme^{-3at/R^{3}}$, which is an exponential distribution falling into the class of Weibull distributions. At $t=40$ (see Fig. 4b) the density of the minimum still qualitatively deviates from an exponential density; while the exponential density is a convex function of $a$ the resulting curve from Fig. 4b for $t=40$ clearly did not yet reach a convex shape in $a$. While the Ornstein-Uhlenbeck process shows a Gumbel distribution in the limit $t\to\infty$, the Bessel process provides an example in which the extreme value distribution falls into the class of Weibull distributions.

4 Concluding perspectives

We used the link between first passage and extremum functionals of reversible ergodic Markov processes in order to formulate the probability density of extreme values in terms of the first passage times. We pushed the connection between these two functionals even further, by utilizing the duality between first passage and relaxation processes [38, 39], which allowed us to determine the statistics of extremes from transition probability densities describing the relaxation towards equilibrium. In their present form our results hold for diffusion in effectively one-dimensional potential landscapes that are sufficiently confining to allow for a discrete eigenspectrum of the corresponding Fokker-Planck operator. Our findings provide a new and deeper perspective on the study of extrema of asymmetric diffusion processes beyond a constant drift. We emphasize that the full probability density of extreme values ($p^{\rm max}$ or $p^{\rm min}$) on arbitrary time-scales still requires the knowledge of the eigenspectrum of the Fokker-Planck operator.

To avoid an eigendecomposition of the Fokker-Planck operator entirely, we established the long time asymptotics of the distribution of extreme values, $p^{\rm max}_{\rm LD}\simeq p^{\rm max}_{\rm LD}$ (or $p^{\rm min}_{\rm LD}\simeq p^{\rm min}$), which accounts for the slowest decaying mode $\propto\rme^{-\mu_{1}t}$ ignoring all faster decaying contributions ($\propto\rme^{-\mu_{2}t},\rme^{-\mu_{3}t},$ etc.). In this large deviation limit we determined explicit bounds on the exact slowest time-scale $\mu_{1}^{-1}$, and showed that these asymptotically tightly bound $\mu_{1}$ from above and from below, which is the central result of this paper.

We illustrated the usefulness of our results by analyzing the statistics of maximum value of the Ornstein-Uhlenbeck process and the minimal distance to the origin of a confined 3d Brownian motion (Bessel process). Our examples underline that the large deviation limit, albeit designed to be asymptotically exact for infinitely long times, approximates the density of the maximum surprisingly well even on relatively short times comparable to the relaxation time, $t\gtrsim 1/\lambda_{1}$. Since $t=1/\lambda_{1}$ reflects the time-scale on which the process tends to decorrelate from the initial condition, the present results describe the statistics of extrema in presence of weak but non-vanishing correlations, and hence go beyond the three classes of limit laws for non-correlated random variables, i.e. the Gumbel, Fréchet, and Weibull distributions [16, 17, 18, 19, 20, 14] as demonstrated on hand of the Ornstein-Uhlenbeck and Bessel process. More generally, it would be interesting to systematically investigate the effect of the potential shape on the limiting extreme value distribution as in Ref. [18].

The remarkable accuracy of the approximation can readily be explained by the interlacing of first passage and relaxation time-scales ($\mu_{1}\leq\lambda_{1}\leq\mu_{2}\ldots$) [38, 39], which renders all higher contributions ($\propto\rme^{-\mu_{2}t},\rme^{-\mu_{3}t}\ldots$) negligibly small compared to the large deviation limit $\propto\rme^{-\mu_{1}t}$ once the condition $t\gtrsim\lambda_{1}^{-1}$ is met.

Our results can be extended and generalized in various ways. Extending the formalism presented here to systems obeying a discrete state Master equation would be straightforward; for example, our main result Eq. (25) would still hold by formally replacing the probability densities from Eq. (24) by the corresponding state probabilities [39]. Interesting and challenging extensions could include the consideration of trapping times [55, 56], spatial disorder [57] and multi-channel transport [58], as well as extreme value statistics in discrete-state Markov processes with a broken time-reversal symmetry [34, 35].

The financial support from the German Research Foundation (DFG) through the Emmy Noether Program “GO 2762/1-1” (to AG) is gratefully acknowledged.