Living on the edge of instability

TL;DR

This paper analyzes the statistical behavior of trajectories in highly unstable potentials, distinguishing between diverging and non-diverging paths, and introduces two limit distributions that generalize Gibbs distribution for such systems.

Contribution

It introduces a novel framework using Q-processes and quasi-stationary distributions to describe non-diverging trajectories in unstable potentials, revealing fundamentally different distribution properties.

Findings

Distribution for diverging trajectories is light-tailed and localized.

Distribution for non-diverging trajectories is heavy-tailed with a power-law decay.

Results apply broadly to nonlinear dynamical models with meta-stable states.

Abstract

Statistical description of stochastic dynamics in highly unstable potentials is strongly affected by properties of divergent trajectories, that quickly leave meta-stable regions of the potential landscape and never return. Using ideas from theory of Q-processes and quasi-stationary distributions, we analyze position statistics of non-diverging trajectories. We discuss two limit distributions which can be considered as (formal) generalizations of the Gibbs canonical distribution to highly unstable systems. Even though the associated effective potentials differ only slightly, properties of the two distributions are fundamentally different for all highly unstable system. The distribution for trajectories conditioned to diverge in an infinitely distant future is localized and light-tailed. The other distribution, describing trajectories surviving in the meta-stable region at the instant of conditioning, is heavy-tailed. The exponent of the corresponding power-law tail is determined by the leading divergent term of the unstable potential. We discuss different equivalent forms of the two distributions and derive properties of the effective statistical force arising in the ensemble of non-diverging trajectories after the Doob h-transform. The obtained explicit results generically apply to non-linear dynamical models with meta-stable states and fast kinetic transitions.

Equations

$$\frac{\mathrm{d}{X}_t}{\mathrm{d}t}=-\frac{1}{\gamma }{V}^{\prime }\left(X_t\right)+\sqrt{2D}{\xi }_t,$$

$$V(x)=\frac{\mu }{3}{x}^3.$$

$$X_t=-{\left[\frac{1}{|y{|}^{2-n}-(n-2)\mu t}\right]}^{\frac{1}{n-2}},n>2.$$

$${\hat{\mathcal{L}}}^{†}=D{\partial }_{xx}^2+{\gamma }^{-1}{\partial }_x{V}^{\prime }(x),$$

$$0<{\lambda }_0<{\lambda }_1<\dots ,$$

$${\hat{\mathcal{L}}}^{†}{p}_n(x)=-{\lambda }_n{p}_n(x).$$

$$\Delta ={\lambda }_1-{\lambda }_0,$$

$${\partial }_{t}P(x,t|y)={\hat{\mathcal{L}}}^{†}P(x,t|y).$$

$$P(x,t|y)=\sum _{n=0}^{\infty }{p}_n(x)s_n(y){\mathrm{e}}^{-{\lambda }_{n}t}.$$

$$\hat{\mathcal{L}}{s}_n(y)=-{\lambda }_n{s}_n(y),$$

$$\hat{\mathcal{L}}=D{\partial }_{yy}^2-{\gamma }^{-1}{V}^{\prime }(y){\partial }_y.$$

$$\int$$

$$S(t|y)=\mathrm{Prob}\{{\tau }_{\mathrm{d}}>t|X_0=y\}=\int P(x,t|y)\mathrm{d}x.$$

$$P(x,t|y)=p_0(x)s_0(y){\mathrm{e}}^{-{\lambda }_{0}t}\left[1+O({\mathrm{e}}^{-\Delta t})\right],$$

$$P(x,t|y)\sim p_0(x)s_0(y){\mathrm{e}}^{-{\lambda }_{0}t},$$

$$S(t|y)\sim s_0(y){\mathrm{e}}^{-{\lambda }_{0}t},$$

$$\hat{\mathcal{H}}$$

$$U(x)=\frac{1}{4D}{\left[\frac{V^{\prime }(x)}{\gamma }\right]}^2-\frac{V^{\prime \prime }(x)}{2\gamma },$$

$$\hat{\mathcal{H}}{\Psi }_n(x)={\lambda }_n{\Psi }_n(x).$$

$$p_0(x)=\frac{1}{\sqrt{Z}}{\mathrm{e}}^{-\beta V(x)/2}{\Psi }_0(x),$$

$$s_0(x)=\sqrt{Z}{\mathrm{e}}^{\beta V(x)/2}{\Psi }_0(x),$$

$$Z$$

$$\int {\mathrm{e}}^{-\beta V(x)}{s}_0(x)\mathrm{d}x=\int {\mathrm{e}}^{-\beta V(x)}{s}_0^2(x)\mathrm{d}x.$$

$$\int$$

$$p_0(x)=\frac{1}{Z}{\mathrm{e}}^{-\beta V(x)}{s}_0(x),$$

$$s_0(x)p_0(x)=\frac{1}{Z}{\mathrm{e}}^{-\beta V(x)}{s}_0^2(x).$$

$$\pi (x,t|y)=\underset{\tau \to \infty }{\lim }{⟨\delta \left(X_t-x\right)⟩}_{{\tau }_{\mathrm{d}}>\tau }=\underset{\tau \to \infty }{\lim }\frac{\int P(z,\tau ;x,t|y)\mathrm{d}z}{S(\tau |y)}.$$

$$\frac{\int P(z,\tau -t|x)P(x,t|y)\mathrm{d}z}{S(\tau |y)}\sim \frac{\int p_0(z)s_0(x){\mathrm{e}}^{-{\lambda }_0(\tau -t)}P(x,t|y)\mathrm{d}z}{s_0(y){\mathrm{e}}^{-{\lambda }_0\tau }},\tau \to \infty .$$

$$\pi (x,t|y)=\frac{s_0(x)}{s_0(y)}{\mathrm{e}}^{{\lambda }_{0}t}P(x,t|y).$$

Full text

Keywords: Unstable dynamics, divergent trajectories, conditional process, quasi-stationary distribution, Q-process, Doob’s h-transform

1 Introduction

Unstable stochastic dynamics is ubiquitous in non-linear models. Its characteristic feature is fast divergence of trajectories that have left a meta-stable state. The divergence restricts precision and duration of experiments [1, 2], and limits applicability of the standard statistical analysis based on the averages, because the latter rapidly diverge with time [3, 4, 5].

In recent theoretical and experimental works [5, 2], the standard description based on moments has been replaced by analysis of local characteristics of the system given by the mode and the curvature near the maximum of the probability density function. The second important point was the shift of attention from the complete ensemble of trajectories to the conditional statistics of non-divergent ones. This allowed to obtained a detailed picture of the unstable dynamics. The theory has been applied to the paradigmatic case of the one-dimensional overdamped Brownian motion

[TABLE]

in the highly unstable cubic potential

[TABLE]

Above, the standard Gaussian white noise satisfies $\langle\xi_{t}\rangle=0$, $\langle\xi_{t}\xi_{t^{\prime}}\rangle=\delta(t-t^{\prime})$, and the diffusion coefficient $D$ is proportional to the thermal energy $k_{\rm B}T$ in accordance with the Einstein relation $D=k_{\rm B}T/\gamma$, where $\gamma$ stands for the friction constant.

The main focus of the aforementioned works has been on the short-time dynamics, zero-noise limit and the so-called quasi-stationary distribution in highly unstable systems. In the present work, we use ideas from the probabilistic theory of conditioned stochastic processes [6, 7] and derive generic behavior of non-diverging trajectories.

First, in Sec. 2, we define precisely what we understand by highly unstable dynamics, demonstrate typical features of trajectories, and emphasize that their divergence is related to the loss of normalization of the propagator. In Sec. 3, we show that highly unstable systems have a discrete spectrum and prove an identity for the generalized partition function. The latter allows us to derive tails of the both discussed limit distributions (Sec. 5.1 and Sec. 6), and the magnitude of an effective statistical force arising in the conditioned ensemble from the Doob h-transform [8] (5.2).

The quasi-stationary distribution (Sec. 6) and the limit distribution of the trajectories that never diverge (Sec. 5) can be regarded as generalizations of the Gibbs canonical distribution to appropriate ensembles of highly unstable systems. Both these distributions reduce to the Gibbs canonical one for stable potentials. Analysis of their properties generic for highly unstable dynamics, particularly based on the notion of effective potentials, are the main topic of the present work.

Alternative methods that can be successfully applied to study highly unstable dynamics include analysis of first-passage times [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 1], and of the so-called nonlinear relaxation times [25, 26]. The analysis of times yields e.g. transition rates for unstable processes. Our present approach is complementary to these methods in the sense that it provides a detailed picture of the particle position statistics in different naturally arising conditional ensembles.

2 Highly unstable dynamics

By highly unstable dynamics we shall understand an overdamped Brownian motion in a potential $V(x)$ decreasing towards $-\infty$ for large $|x|$ so fast that the particle reaches (minus) infinity in a finite time. We shall call such potentials as highly unstable. Equivalently, they can be defined as the class of potentials for which the mean first-passage time from any point $x$ to $x=\pm\infty$ is finite [24].

The simplest examples of highly unstable potentials include the cubic, and the inverted quartic potentials. Generally speaking, highly unstable potentials decrease towards minus infinity at least as $-|x|^{n}$ for large $|x|$, with $n>2$. To see that such potentials are highly unstable, we inspect the zero-noise limit of the Langevin equation (1) with $V(x)\sim-|x|^{n}$ when $x\to-\infty$. For the initial condition $X_{0}=y$ placed into the unstable region, $y<0$, we obtain the deterministic solution

[TABLE]

This solution reaches minus infinity at the finite time $t_{\infty}=|y|^{2-n}/[(n-2)\mu]$, when the denominator attains zero. Notice that for the inverted parabolic potential ($n=2$), we get the exponential divergence of $X_{t}$. Even though the exponential divergence can be considered as “fast” in some contexts, it is much “slower” than the divergence generated by Eq. (3). Therefore, the inverted parabolic potential is not highly unstable.

Equivalently, one can determine the condition $n>2$ starting from the exact expression for the mean-first passage time to (minus) infinity. The mean time is known analytically in one dimension [27]. For the detailed discussion related to the high instability of the cubic potential see e.g. Ref. [24] and references therein.

Figure 1(b) shows a typical sample of trajectories starting at $y=0$ and diffusing in the highly unstable cubic potential. A trajectory in the sample can consist of three distinct segments. First, on the potential plateau, i.e., in the region near the origin approximately determined by the condition that the cubic potential is weak compared to thermal energy, $|V(x)|/k_{\rm B}T<1$, the trajectory undergoes the (almost) free diffusion. Second, occasionally, the trajectory is reflected by the strong cubic repulsion and returns quickly to the plateau if it reaches a large positive $x$ [upward spikes in Fig. 1(b)]. Third, the last part of the trajectory begins if the particle moves left from the plateau and reaches highly unstable region $x<0$. There, the decreasing branch of the cubic potential is extremely strong compared to the thermal noise, $|V(x)|/k_{\rm B}T\gg 1$, and the trajectory is quickly dragged towards minus infinity. The divergence is nearly deterministic, in accord with Eq. (3) for $n=3$.

A characteristic feature of highly unstable dynamics, which we prove in Sec. 3, is that the Fokker-Planck operator corresponding to the Langevin equation (1),

[TABLE]

has a discrete spectrum of eigenvalues

[TABLE]

[TABLE]

That is, for highly unstable dynamics, there exists an isolated leading decay rate $\lambda_{0}$, and the strictly positive gap

[TABLE]

between the eigenvalues corresponding to the slowest and the second-slowest decaying modes.

The propagator, or the probability density function (PDF) of the particle position at time $t$, given that initially the particle was at $y$, satisfies the Fokker-Planck equation

[TABLE]

The eigenvalues determine decay rates of individual modes with amplitudes given by the eigenfunctions $p_{n}(x)$ defined in Eq. (6). The propagator can be expressed in terms of the eigenfunctions as [28]

[TABLE]

Above, the functions $s_{n}(y)$ depend on the initial particle position, $y$, and follow from the adjoint eigenvalue problem

[TABLE]

with the operator

[TABLE]

adjoint to the Fokker-Planck operator $\hat{\mathcal{L}}^{\dagger}$. The adjoint operator is called the generator of the Markov process. The two sets of eigenfunctions satisfy the following normalization conditions

[TABLE]

We shall refer to $p_{n}(x)$ and $s_{n}(y)$ as the right and left eigenfunctions of the Fokker-Planck operator $\hat{\mathcal{L}}^{\dagger}$, respectively.

The high instability of dynamics manifests itself through exponential decay of normalization $S(t|y)$ of the propagator $P(x,t|y)$ (9). The normalization equals to the probability that a trajectory has not diverged up to time $t$ and thus we call it as the survival probability [9]. In symbols, we have

[TABLE]

Above, $\tau_{d}$ stands for the time of divergence of a trajectory. According to Eq. (9), in the long-time limit only the most-slowly decaying term contributes significantly to the sum, i.e.,

[TABLE]

which we will simply write as

[TABLE]

and

[TABLE]

where “$\sim$” means that the ratio of the two expressions converges to unity with increasing time.

Physical meaning of the exponential decay (17) can be understood based on the illustration in Fig. 1(b). The decay arises from the fast divergence of majority of the trajectories which reach $X_{t}=-\infty$ in a finite time. They are not included in the statistics described by the propagator $P(x,t|y)$ for any finite $x$, $x\in(-\infty,+\infty)$, and their weight is given by $[1-S(t|y)]$.

Finally, notice that for a confining potential $V_{\rm c}(x)$, the dynamics becomes stable and conservative, and we have $\lambda_{0}=0$, $s_{0}(x)=1$, and $s_{0}(x)p_{0}(x)=p_{0}(x)=\pi_{\rm eq}(x)={\rm e}^{-\beta V_{\rm c}(x)}/Z$, where the equilibrium partition function ensures the normalization condition, $Z=\int{\rm e}^{-\beta V_{\rm c}(x)}{\rm d}x$. Contrary to this, in the highly unstable non-conservative case, the eigenfunction $s_{0}(x)$ is no-longer constant. Below, we shall see that the expressions $Q_{\rm st}(x)=p_{0}(x)$ and $\pi_{\rm st}(x)=s_{0}(x)p_{0}(x)$, which are identical in the stable case, play roles of stationary distributions in proper conditioned ensembles of non-diverging trajectories. The two ensembles, described by these distributions, have fundamentally different statistical properties, cf. Sec. 5.1 and Sec. 6.

3 Schrödinger equation: Normalization of eigenvectors

The eigenproblem for the non-Hermitian Fokker-Planck operator involves distinct left and right eigenfunctions for each eigenvalue $-\lambda_{n}$. The two sets of eigenfunctions satisfy different boundary conditions and should be normalized properly to form orthogonal spaces, as we have described above in Sec. 2.

The non-Hermitian eigenproblem can be transformed into the regular Hermitian one by a similarity transformation that recasts the Fokker-Planck operator (or equivalently the generator) into the Hermitian Hamiltonian [29, 27, 30]:

[TABLE]

where $\beta=1/k_{\rm B}T$, and the diffusion coefficient, $D$, determines the inverse mass of the representative quantum particle. The transformed potential,

[TABLE]

is confining for all highly unstable potentials $V(x)$. For the cubic potential, $U(x)$ is a quartic well [the term proportional to $(V^{\prime}(x))^{2}$] with a slight asymmetry and shift of the minimum caused by the linear term arising from $V^{\prime\prime}(x)$, see Fig. 2.

The transformation provides two valuable insights. The first results from the discrete structure of the spectrum of $\hat{\mathcal{H}}$ that is equivalent to the discrete well-separated spectrum of $\hat{\mathcal{L}}^{\dagger}$. Indeed, the spectrum of $\hat{\mathcal{H}}$ is identical to that of $\hat{\mathcal{L}}^{\dagger}$ up to the sign:

[TABLE]

Because $U(x)$ is always confining for highly unstable $V(x)$, the spectrum is always discrete as it is the case for all confining infinite potential wells [31]. Therefore, the inequalities (5) follow directly from the high instability of $V(x)$.

The second insight arises from the transformation between the wave functions $\Psi_{n}(x)$ and eigenfunctions of $\hat{\mathcal{L}}^{\dagger}$. The ground state wave function, associated with the energy $\lambda_{0}$, is related to the eigenfunctions $p_{0}(x)$ and $s_{0}(x)$ of $\hat{\mathcal{L}}^{\dagger}$ via

[TABLE]

respectively. The factor $Z$ is independent of $x$ and it plays no role in the construction of the propagator (9) for the diffusion problem. Hence, $Z$ is usually omitted in standard literature, e.g. [29, 27, 30], when discussing the spectral solution of the Fokker-Planck equation. On the contrary, the factor $Z$, to which we shall refer as to the generalized partition function, is fundamental for the understanding of the Q-process.

Two equivalent integral representations of the partition function, used below to explore properties of stationary distributions, are

[TABLE]

Hence, we obtain the identity

[TABLE]

Even though the functions ${\rm e}^{-\beta V(x)}s_{0}(x)$ and ${\rm e}^{-\beta V(x)}s_{0}^{2}(x)$ are strikingly different, they enclose equal areas. Figure 2 shows both the localized function ${\rm e}^{-\beta V(x)}s_{0}^{2}(x)$ with light tails, and the function ${\rm e}^{-\beta V(x)}s_{0}(x)$ possessing a heavy tail for negative $x$.

The identity (18wx) follows from the normalization conditions (12) and (13), which, for the eigenfunctions $p_{0}(x)$ and $s_{0}(x)$, read

[TABLE]

The combination of Eqs. (18u) and (18v) yields

[TABLE]

The requirement that right-hand sides of both equations are normalized to unity immediately gives us the identity (18wx). In other words, the constant $Z$ in Eqs. (18u) and (18v) is crucial for consistency of Eqs. (18waa) and (18wab), relating the eigenfunctions $p_{0}(x)$ and $s_{0}(x)$, with the normalization conditions (18wy) and (18wz).

For the Gibbs canonical equilibrium distribution in a confining potential $V_{\rm c}(x)$, the generalized partition function $Z$ is identical to the equilibrium partition function. In equilibrium, we have $\lambda_{0}=0$, $s_{0}(x)=1$, and the corresponding right eigenfunction of $\hat{\mathcal{L}}^{\dagger}$ is the Gibbs canonical distribution $p_{0}(x)={\rm e}^{-\beta V_{\rm c}(x)}/Z$, with $Z=\int{\rm e}^{-\beta V_{\rm c}(x)}{\rm d}x$.

4 Q-process: Dynamics

Figure 1(c) shows three distinct dynamical regimes observed on trajectories that do not diverge in a highly unstable potential $V(x)$. In the transient regime I, the relaxation of the initial condition takes place. In the stationary regime II, the initial condition is already forgotten and the ensemble of trajectories resembles sample paths of a stationary conservative process diffusing away from the instability region. Within II, the particle position is described by the time-independent PDF $\pi_{\rm st}(x)$. Finally, in the transient regime III, the distribution departs from $\pi_{\rm st}(x)$ and approaches the heavy-tailed quasi-stationary distribution $Q_{\rm st}(x)$ discussed in Sec. 6.

Durations of the transient regimes I and III can be estimated by the inverse difference of decay rates $1/\Delta$, $\Delta=\lambda_{1}-\lambda_{0}$, that is finite due to high instability of $V(x)$. The finite transient times allow for an exact statistical description of trajectories in the regimes I and II. Formally, this is done by shifting the last transient regime III towards infinite times. This limit has no significant influence on trajectory statistics within regimes I and II. Remarkably, the limit transforms the non-Markovian process, conditioned on a non-local property (non-divergence up to a given time), into a Markovian one, that can be analyzed using standard tools as the Langevin and the Fokker-Planck equation.

Let us choose an “intermediate” time $t$, and the “final” time of conditioning $\tau$, such that $0\leq t\leq\tau$. In Fig. 1(c), we have $\tau=15$ and $t\in[0,15]$. All trajectories of interest are assumed to diverge after time $\tau$. That is, the time of divergence $\tau_{\rm d}$ of any trajectory satisfies $\tau_{\rm d}>\tau$. Now, we can formally shift the conditioning time $\tau$ to infinity and define the distribution of trajectories at finite time $t$:

[TABLE]

The limit $\tau\to\infty$ means that the process is conditioned to never diverge. In mathematical literature, such a process is commonly denoted as the Q-process [32, 33].

The propagator $\pi(x,t|y)$ for the Q-process is derived by a direct evaluation of the above limit. Employing the Markov property, we rewrite the joint probability density in (18wac) as the product of propagators for the intervals $[0,t]$ and $[t,\tau]$: $P(z,\tau;x,t|y)=P(z,\tau-t|x)P(x,t|y)$. Then, for large $\tau$ and fixed $t$, we get

[TABLE]

The conditioning time $\tau$ cancels out and the integral over $z$ equals to one. The limit shifts the time of conditioning and the whole transient region III in Fig. 1(c) towards infinite times and removes information on III from the dynamics of trajectories at any fixed $t$.

The right-hand side of Eq. (18wad) yields explicit form of the propagator for the Q-process:

[TABLE]

It is straightforward to verify that the Q-process is a conservative Markov process. The normalization condition, $\int\pi(x,t|y){\rm d}x=1$, is verified most easily using the eigenfunction expansion of $P(x,t|y)$ (9) and the orthogonality relation for the eigenfunctions (13). The Chapman-Kolmogorov equation, $\pi(x,t_{1}+t_{2}|y)=\int\pi(x,t_{2}|z)\pi(z,t_{1}|y){\rm d}z$, follows directly from the Chapmann-Kolmogorov equation for the bare propagator $P(x,t|y)$.

The Q-process is the driven diffusion process and, as such, its propagator $\pi(x,t|y)$ satisfies the backward equation

[TABLE]

with the generator

[TABLE]

The generator (18waga) is obtained directly by taking the time derivative of Eq. (18wae), and using that $s_{0}^{-1}(\hat{\mathcal{L}}\,s_{0})=-\lambda_{0}$. The explicit form (18wagb) follows from (18waga) after some algebra. The Fokker-Planck operator adjoint to the generator $\hat{\mathcal{L}}_{\rm s}$ reads

[TABLE]

Equivalently, the dynamics of trajectories can be described by the Langevin equation corresponding to the Fokker-Planck operator (18wagahb),

[TABLE]

In mathematics, the similarity transformation (18waga) is known as the Doob h-transform [8]. It leaves intact the diffusion coefficient and introduces an effective drift that keeps the trajectories away from the unstable region of the potential. For a more rigorous discussion of the Q-process, we refer to the extensive mathematical literature on the topic, see e.g. [7] and references therein. Here, our aim was to introduce these concepts heuristically and accessibly for a wide audience of physicists.

In physics, an analogue of the Q-process can be found in theory of large deviations as one of effective processes that can be used to generate ensembles associated with rare values of time-extensive observables [34, 35, 36, 37, 38, 39]. For a comprehensive review of Markov processes conditioned on large deviations see Ref. [6], where also the connection to the original Doob’s conditioning is described in a great detail.

5 Q-process: Limit distribution

For any initial condition, the dynamics of Q-process converges towards the stationary regime II with the relaxation time $1/\Delta$. The limit distribution of the Q-process can be computed as the long-time limit of its propagator,

[TABLE]

The first equality shows that any trace of the initial condition $y$ is forgotten. The second one emphasizes importance of correct ordering of the two limits and suggests that the limit distribution $\pi_{\rm st}(x)$ differs from the quasi-stationary distribution $Q_{\rm st}(x)$. The latter appears at the final instant $\tau$ of the transient regime III, which we have shifted to infinity before taking the limit $t\to\infty$, cf. Eq. (18wac) and Eq. (18wagahbe) below.

The limit distribution $\pi_{\rm st}(x)$ describes the stationary state of trajectories which will diverge in (an infinitely) distant future. The explicit evaluation of the $t\to\infty$ limit with the aid of the leading asymptotic behavior of the propagator (16) yields

[TABLE]

Hence the limit distribution is given by

[TABLE]

The PDF $\pi_{\rm st}(x)$ is the stationary distribution (corresponding to the eigenvalue 0) of the Fokker-Planck operator (18wagahb). The equality $\hat{\mathcal{L}}^{\dagger}_{\rm s}\pi_{\rm st}=0$ can be verified directly by applying the Fokker-Planck operator (18wagaha) on the PDF (18wagahal).

The relations (18u) and (18v), between the eigenfunctions $p_{0}(x)$, $s_{0}(x)$ and the ground state wave function $\Psi_{0}(x)$ reveal the physical meaning of $\Psi_{0}(x)$:

[TABLE]

The equality implies that powerful quantum-mechanical methods for approximating ground state wave functions can be applied to discuss $\pi_{\rm st}(x)$ for highly unstable dynamics.

Another expression for $\pi_{\rm st}(x)$ follows from the relation (18wab) between the eigenfunctions $p_{0}(x)$ and $s_{0}(x)$ established in Sec. 3,

[TABLE]

In this form, $\pi_{\rm st}(x)$ resembles the Gibbs canonical distribution with the density of states $g(V){\rm d}V=s_{0}^{2}(x){\rm d}x$, and $Z=\int{\rm e}^{-\beta V}g(V){\rm d}V$. The divergence of the Boltzmann factor ${\rm e}^{-\beta V}$ for $x\to-\infty$ is compensated by the fast decrease of the squared eigenfunction $s_{0}^{2}(x)$, accounting for the fact that the non-diverging trajectories do not populate the unstable region $x\ll 0$, and thus the PDF $\pi_{\rm st}(x)$ is normalizable.

Last but not least, it is convenient to reshape Eq. (18wagahan) into the Gibbs canonical form containing an effective potential $V_{\rm eff}$. The result is

[TABLE]

where

[TABLE]

The minus derivative of $V_{\rm s}(x)$ gives us the effective force,

[TABLE]

arising from the tendency of non-diverging trajectories to avoid the highly unstable region $x\ll 0$. In the Langevin equation (18wagahai) and in the corresponding Fokker-Planck operator (18wagahb), the force (18wagahaq) acts in a superposition with the actual external drift force given by $-V^{\prime}(x)$.

Even though the limit distribution of the non-diverging trajectories (18wagahao) and the generalized partition function (18wb) are formally similar to their counterparts in the equilibrium canonical ensemble, there is an important physical difference. In the equilibrium case, both $\pi_{\rm st}(x)$ and $Z$ are determined solely by the temperature and the external potential $V(x)$. On the contrary, $\pi_{\rm st}(x)$ and $Z$ for the Q-process depend also on the friction constant $\gamma$ through the eigenfunction $s_{0}(x)$.

Hence, $\pi_{\rm st}(x)$ and $Z$ for non-diverging trajectories are kinetic quantities whose values depend on time scales involved in the problem [and not only on the temperature and the external potential $V(x)$]. Qualitatively similar situation appears in theory of non-equilibrium steady states of processes driven by e.g. couplings to several reservoirs with different chemical potentials and temperatures, or/and by non-conservative external forces. In such steady states, kinetic parameters like friction coefficients (called also non-dissipative or frenetic) often determine the overall shape of the stationary probability distribution and magnitudes of the corresponding currents [40, 41, 42, 43, 44].

Analogous strong dependencies on kinetic parameters can be expected for the limit distribution of non-diverging trajectories in highly unstable dynamics. The two cases exhibit an interesting difference in the mathematical structure of the probability distributions. In the non-equilibrium steady states, we have $s_{0}(x)=1$ and both the kinetic and the thermodynamic effects are mixed in a single PDF given by $p_{0}(x)$. On the other hand, in the highly unstable dynamics, the two effects are separated and contained within the two contributions to the effective potential (18wagahap). Further investigation in this direction, focusing in particular on higher-dimensional models, may reveal new interesting phenomena for conditioned processes.

5.1 Tails of the limit distribution and effective potential

For large $x$ outside the unstable region, the eigenfunction $s_{0}(x)$ [and the potential $V_{\rm s}(x)$] approaches a constant and the force $F_{\rm s}(x)$, repelling non-diverging trajectories from the highly unstable region, vanishes. Hence the right tail of the limit distribution $\pi_{\rm st}(x)$ can be approximated by the Boltzmann tail:

[TABLE]

In this sense, the force $F_{\rm s}(x)$ is short-ranged and has no significant effect on dynamics for $x\gg 0$, where the effective potential in Eq. (18wagahao) is given by the actual external potential, $V_{\rm eff}(x)\sim V(x)$, for $x\gg 0$.

The left tail of the limit distribution can be inferred directly from the definition of $Z$ in Eq. (18wa) and the eigenvalue equation for $s_{0}(x)$,

[TABLE]

We multiply the equation with the Boltzmann factor ${\rm e}^{-\beta V}$ and integrate over $x$. Equation (18wa) assures that the integral is convergent and $Z$ appears on the left-hand side. After the integration, we obtain

[TABLE]

Therefore, the controlling asymptotic factor of $s_{0}(x)$ for $x\ll 0$ is the exponential ${\rm e}^{\beta V}$. It is straightforward to verify that the above limit is satisfied by the function

[TABLE]

The limit distribution thus decays exponentially as

[TABLE]

the exponential decay being modified by a polynomial prefactor determined by the unstable part of $V(x)$, $x\to-\infty$. Note that the constant prefactor in this equation is exact because it follows from the limit on the right-hand side of Eq. (18wagahat).

The result (18wagahav) provides complete information on the asymptotic behavior of the total effective potential $V_{\rm eff}(x)=V(x)+V_{\rm s}(x)$ in Eq. (18wagahao) and in particular on its part $V_{\rm s}(x)$ arising due to the conditioning. The latter behaves as

[TABLE]

The first term on the right-hand side compensates the external unstable potential $V(x)$ and contributes by the exponent $\beta V(x)$ responsible for the exponential decay in Eq. (18wagahau). The second term creates a logarithmic barrier $-2k_{\rm B}T(n-1)\log(x)$. Interestingly enough, the barrier is impenetrable for all cases of interest in the present paper where $n>2$; for a further discussion see classification of boundary conditions for the Bessel process, e.g. in Refs. [45, 46]. Finally, we emphasize that the relation (18wagahaw) holds only asymptotically, i.e., the effective potential has no singularity around $x=0$ where majority of non-diverging trajectories live.

5.2 Mean squared force

Having understood the asymptotic behaviors of the effective potential $V_{\rm eff}(x)$, it is interesting to ask, what is an average magnitude of the force $F_{\rm s}(x)$ and how it depends on the model parameters $T$, $\gamma$, $\mu$, and $n$. A simple answer is provided in this Section. It is based on the second expression for $Z$ in Eq. (18wb) and on a scaling analysis.

Intuitively we expect the magnitude of the force $F_{\rm s}(x)$ to increase with the degree of instability. Indeed, it turns out that the mean squared amplitude of the force is proportional to the decay rate

[TABLE]

The average is taken with respect to $\pi_{\rm st}(x)$. Equation (18wagahax) follows directly from the differential equation (18wagahas) for $s_{0}(x)$, and from Eq. (18wa) for the partition function containing the second power of $s_{0}(x)$. We multiply both sides of Eq. (18wagahas) by $s_{0}(x){\rm e}^{-\beta V(x)}$, integrate over $x$, identify $Z$ on the right-hand side in accordance with Eq. (18wb), and get

[TABLE]

where we have performed the per partes integration on the right-hand side. After dividing Eq. (18wagahay) with $Z$ and adjusting properly the $k_{\rm B}T$ factor, we get Eq. (18wagahax).

Equation (18wagahax) is valid generally for any potential $V(x)$. Its further explicit analysis is possible for monomial unstable potentials $V(x)=\mu x^{n}/n$ via the scaling analysis. Defining the time and the length scales as [2]

[TABLE]

respectively, the eigenvalue problem (6) attains the dimensionless form

[TABLE]

where $\zeta=x/\tilde{x}$ stands for the transformed coordinate variable, and the transformed dimensionless decay rate $\theta_{0}=\lambda_{0}\tilde{t}$, depends on the order $n$ only. Similar transformation can be applied to the adjoint eigenvalue problem for the left eigenfunctions.

Rescaling the position variable by $\tilde{x}$ and the decay rate with $\tilde{t}$, we obtain from the exact equation (18wagahax) the explicit dependence

[TABLE]

of the force magnitude on the model parameters. Similar dependence holds for the mean value of the force,

[TABLE]

where the prefactor $A$ depends on $n$ only.

The effective force always increases with the temperature $T$ and with the potential amplitude $\mu$. The dependence on $\mu$ becomes weaker for larger powers $n$, which suppress the potential magnitude at the plateau around $x=0$, where majority of non-diverging trajectories live. For large $n$, the temperature dependence in (18wagahbc) becomes linear for the same reason.

6 Quasi-stationary distribution

In this section, we compare the limit distribution of the Q-process $\pi_{\rm st}(x)$ and the quasi-stationary distribution $Q_{\rm st}(x)$. The latter arises at the end of the transient regime III in Fig. 1(c). Specifically, we shall use the derived properties of the effective potential and relations between the eigenvectors to illustrate the most striking difference between the two distributions.

The statistics of trajectories at time $\tau$ that are non-diverging at least up to time $\tau$ is characterized by the conditional PDF

[TABLE]

where the inequality $\tau_{\rm d}>\tau$ means that each trajectory in the conditioned ensemble diverges later than at $\tau$. In contrast to Eq. (18wac), here, the time of conditioning $\tau$ is identical with the time of observation of the process. Existence of the spectral gap $\Delta$ implies that the conditioned PDF (18wagahbd) converges exponentially fast to the limit

[TABLE]

The limit is known as the quasi-stationary distribution in mathematical literature, where it has first been introduced in context of branching processes and absorbed Marov chains [47, 48, 49, 50] and since then has been studied extensively in context of population dynamics with extinction [7]. The extinction event in population dynamics is in our context represented by the trajectory divergence.

Evaluation of the limit (18wagahbe) with the aid of the long-time expressions  (16) and (17) for $P(x,\tau|y)$ and $S(\tau|y)$, respectively, shows that it is equal to the eigenfunction $p_{0}(x)$, i.e.,

[TABLE]

Then, Eq. (18waa) can be used to recast this relation into the form

[TABLE]

analogous to the one we have examined for $\pi_{\rm st}(x)$, cf. Eq. (18wagahan). The two equations differ just by a power of $s_{0}(x)$ on their right-hand sides. The generalized partition functions $Z$ are equal in the both equations because of the identity (18wx), expressing that areas enclosed by nominators in Eqs. (eq:Qsts0) and (18wagahan) are the same even though the two functions are rather different, as illustrated in the right panel of Fig. 2.

To pursue the comparison further, the quasi-stationary distribution (18wagahbg) can be written in the form

[TABLE]

resembling the canonical-like expression for $\pi_{\rm st}(x)$ in Eq. (18wagahao). There is just a minor difference between expressions for the two distributions given by the factor $1/2$ in the above exponent, i.e., $Q_{\rm st}(x)={\rm e}^{\beta V_{\rm s}(x)/2}\pi_{\rm st}(x)$. This minor difference has no effect on the right tail of the quasi-stationary distribution. Hence the tail is the same as that of $\pi_{\rm st}(x)$ in Eq. (18wagahar),

[TABLE]

because of converges of $s_{0}(x)$ to a constant.

Contrary to the similarity of right tails, the factor $1/2$ in Eq. (18wagahbh) strongly influences the asymptotic behavior for $x\to-\infty$. Based on Eq. (18wagahau), the left tail assumes the form

[TABLE]

Specifically, for $V(x)=\mu x^{n}/n$ the tail decays as the power law

[TABLE]

Hence, in contrast to $V_{\rm eff}(x)=V(x)+V_{\rm s}(x)$ in the Q-process limit distribution (18wagahao), the modified effective potential in the present case of the quasi-stationary distribution (18wagahbh), $[V(x)+V_{\rm s}(x)/2]$, is not strong enough to confine trajectories to the stable region. The power-law decay (18wagahbk), arising from the asymptotically logarithmic potential in Eq. (18wagahbh),

[TABLE]

allows trajectories to populate any negative $x$ exploring the whole unstable region with a slowly decaying probability.

Note that for the cubic potential, no integer moments of the quasi-stationary distribution exist and hence the standard statistical analysis of particle position based on averages provides no information on the actual particle dynamics. An alternative approach based on local quantities (position of the PDF maximum and curvature near the maximum) has been proposed [3, 5] and its utility verified experimentally in [2] on colloidal particles diffusing in a nonlinear optical trap illustrated in Fig. 1(a).

We also note that it is possible to construct an effective conservative Markov process that converges towards the distribution $Q_{\rm st}(x)$ in the long-time limit [5]. The process is based on an appropriate return (or resetting) of the diverging trajectories. The effective process should not be confused with the Q-process discussed above. The two has different limiting distributions and the Q-process describes faithfully the conditioned dynamics whereas the effective process of Ref. [5] converges towards $Q_{\rm st}(x)$ only asymptotically in the long-time limit and its transient dynamics may differ from the actual one.

7 Conclusions and perspectives

We have discussed properties of the position statistics of non-diverging stochastic trajectories diffusing in highly unstable potentials. Two stationary distributions arising in this case can be formally expressed in the Gibbs canonical form with the effective potentials depending on the left eigenvector of the Fokker-Planck operator.

The limit distribution of the Q-process describes stationary statistics of trajectories conditioned to be non-diverging for an infinitely long time. For all highly unstable systems, the distribution is localized with exponentially decaying tails. The effective force arising from the conditioning keeps trajectories away from the unstable region of the potential and its strength increases with the degree of instability.

On the other hand, the quasi-stationary distribution is heavy-tailed towards the instability and hence it allows trajectories to be located anywhere in the unstable region. The slow power-law decay of the quasi-stationary distribution is determined by the derivative of the external potential, cf. Eq. (18wagahbk). Here, an interesting question arises whether the theory of non-normalizable distributions [51] may contribute to our understanding of unstable dynamics with non-integer exponents $n$ smaller than 3.

As shown in Sec. 3, the basic property of highly unstable potentials is that the spectrum of the Fokker-Planck operator is negative and discrete. Other systems satisfying the same property include confined diffusions with absorbing boundaries. At the same time, our results imply that the magnitude of the effective force, which confines conditioned trajectories of the Q-process, decreases with the decreasing instability of the potential, cf. Sec. 5.2. Therefore, it may be interesting to explore the dynamics for other unstable potentials and identify the precise conditions, when the Q-process becomes transient.

Finally, we remark that the similarity of the two distributions with the Gibbs canonical ensemble is purely formal. In general, the distributions for the highly unstable systems, depended on kinetic parameters such as the friction coefficients. This dependence implies a remarkable possibility to control the overall shape of the two distributions by varying the kinetic parameters only. This effect should become more pronounced in higher-dimensional systems and it represents another interesting area for a further investigation.

We gratefully acknowledges financial support by the Czech Science Foundation (project No. 17-06716S). VH also thanks for support by the Humboldt foundation. AR is grateful to organizers of the nesmcq18 conference, where a part of the present work has been discussed. AR acknowledges financial support from the Portuguese Foundation for Science and Technology (FCT) under Contracts nos. PTDC/FIS-MAC/28146/2017 (LISBOA-01-0145-FEDER-028146) and UID/FIS/00618/2019.