Reaction rates and the noisy saddle-node bifurcation: Renormalization group for barrier crossing

TL;DR

This paper develops a renormalization group framework to understand barrier crossing times in noisy systems, especially near the saddle-node bifurcation where the barrier vanishes, revealing universal scaling laws and distribution approximations.

Contribution

It introduces a novel renormalization group approach to analyze the noisy saddle-node transition and derives universal scaling behavior for escape times.

Findings

Derived universal scaling laws for mean escape time near the bifurcation.

Developed an approximation for the distribution of barrier escape times.

Linked barrier crossing phenomena with bifurcation and renormalization group theories.

Abstract

Barrier crossing calculations in chemical reaction-rate theory typically assume that the barrier is large compared to the temperature. When the barrier vanishes, however, there is a qualitative change in behavior. Instead of crossing a barrier, particles slide down a sloping potential. We formulate a renormalization group description of this noisy saddle-node transition. We derive the universal scaling behavior and corrections to scaling for the mean escape time in overdamped systems with arbitrary barrier height. We also develop an accurate approximation to the full distribution of barrier escape times by approximating the eigenvalues of the Fokker-Plank operator as equally spaced. This lets us derive a family of distributions that captures the barrier crossing times for arbitrary barrier height. Our critical theory draws links between barrier crossing in chemistry, the renormalization group, and bifurcation theory.

Figures

Equations

$$\dot{x}=f(x)+g(x)\xi (t).$$

$$\frac{dx}{dt}=\sum _{n=0}^{\infty }{ϵ}_n{x}^n+\xi (t)\sum _{n=0}^{\infty }{g}_n{x}^n.$$

$$\frac{d\hat{x}}{d\hat{t}}=\sum _{n=0}^{\infty }{ϵ}_n{b}^{2-n}{\hat{x}}^n+\hat{\xi }(\hat{t})\sum _{n=0}^{\infty }{b}^{3/2-n}{g}_n{\hat{x}}^n.$$

$$\frac{d{ϵ}_n}{d\ell }=(2-n){ϵ}_n,\frac{d{g}_n}{d\ell }=(3/2-n)g_n.$$

$$\tau =g_0^{-2/3}\mathcal{T}\left(\{{ϵ}_n/g_0^{2(2-n)/3}\},\{g_n/g_0^{1-2n/3}\}\right),$$

$$x={\int }^{\tilde{x}}\frac{g_0}{g(y)}dy,$$

$$\tau =g_0^{-2/3}\mathcal{T}({ϵ}_0/g_0^{4/3}).$$

$$\mathcal{T}(\alpha )\sim \frac{\pi }{|\alpha {|}^{1/2}}{e}^{\frac{8}{3}|\alpha {|}^{3/2}},\alpha \to -\infty .$$

$$\mathcal{T}(\alpha )\sim \frac{\pi }{{\alpha }^{1/2}},\alpha \to \infty .$$

$$\mathcal{P}(t)={\int }_{-\infty }^{x_f}\mathcal{P}(x,t)dx.$$

$$\tau (x_i|x_f)={\int }_0^{\infty }tp(t)dt={\int }_0^{\infty }\mathcal{P}(t)dt,$$

$$-\frac{d\mathcal{P}(x,t)}{d{t}_i}=-V^{\prime }(x_i)\frac{d\mathcal{P}(x,t)}{d{x}_i}+\frac{1}{2}{g}_0^2\frac{d^2\mathcal{P}(x,t)}{d{x}_i^2}.$$

$$\frac{1}{2}{g}_0^2{\tau }^{\prime \prime }(x_i|x_f)-V^{\prime }(x_i){\tau }^{\prime }(x_i|x_f)=-1.$$

$$\tau (x_i|x_f)=\frac{2}{g_0^2}{\int }_{x_i}^{x_f}dy{\int }_{-\infty }^{y}dz{e}^{-\frac{2}{g_0^2}\left[V(z)-V(y)\right]},$$

$$\mathcal{T}(\alpha )=2^{1/3}{\pi }^2\left[{\text{Ai}}^2(-2^{2/3}\alpha )+{\text{Bi}}^2(-2^{2/3}\alpha )\right],$$

$$\mathcal{T}(\alpha ,{\chi }_{-},{\chi }_{+})=\mathcal{T}(\alpha )-{\mathcal{T}}_{-}(\alpha ,{\chi }_{-})-{\mathcal{T}}_{+}(\alpha ,{\chi }_{+}).$$

$${\mathcal{T}}_{-}(\alpha ,{\chi }_{-})=2{\int }_{-\infty }^{{\chi }_{-}}dy{\int }_{-\infty }^{y}dz{e}^{-2\left(y^3/3+\alpha y-z^3/3-\alpha z\right)}$$

$${\mathcal{T}}_{+}(\alpha ,{\chi }_{+})=2{\int }_{{\chi }_{+}}^{\infty }dy{\int }_{-\infty }^{y}dz{e}^{-2\left(y^3/3+\alpha y-z^3/3-\alpha z\right)}.$$

$$\begin{array}{rl}{\mathcal{T}}_{\pm }(\alpha ,{\chi }_{\pm })=& |\alpha {|}^{-1/2}{\tan }^{-1}(|\alpha {|}^{1/2}/|{\chi }_{\pm }|)\\ & +\frac{1}{{\chi }_{\pm }^4}\left(\frac{1}{4}-\frac{\alpha }{2{\chi }_{\pm }^2}+\frac{3{\alpha }^2}{4{\chi }_{\pm }^4}\right)\\ & \pm \frac{1}{{\chi }_{\pm }^7}\left(\frac{5}{14}-\frac{13\alpha }{9{\chi }_{\pm }^2}\right)+O({\chi }_{\pm }^{-10}).\end{array}$$

$${\mathcal{T}}_3(\alpha )\stackrel{\alpha \ll 0}{\to }\pi \sqrt{|\alpha |}{e}^{\frac{8}{3}|\alpha {|}^{3/2}}(8|\alpha {|}^{3/2}+11)/8.$$

$$\mathcal{T}(\alpha ,\beta ,\gamma )$$

$$p(t)=\mathcal{E}({\lambda }_0)\ast \mathcal{E}({\lambda }_1)\ast \cdots \ast \mathcal{E}({\lambda }_N),$$

$$\log \tilde{p}(\omega )=\sum _{m=1}^{\infty }{\kappa }_m(i\omega {)}^m/m!$$

$${\kappa }_m=(m-1)!\sum _{n=0}^N\frac{1}{{\lambda }_n^m}.$$

$$\frac{\partial \mathcal{P}(x,t)}{\partial t}=\frac{\partial }{\partial x}\left(V^{\prime }(x)\mathcal{P}(x,t)\right)+\frac{1}{2}{g}_0^2\frac{{\partial }^2\mathcal{P}(x,t)}{\partial x}.$$

$$\frac{\partial\sigma(x,t)}{\partial t}=\left\{\frac{g_{0}^{2}}{2}\frac{\partial^{2}}{\partial x^{2}}+\frac{1}{2}V^{\prime\prime}(x)-\frac{1}{2g_{0}^{2}}\Big{(}V^{\prime}(x)\Big{)}^{2}\right\}\sigma(x,t),$$

$${\lambda }_n=g_0^{2/3}{\Lambda }_n(\alpha ),$$

$${\mathcal{K}}_m(\alpha )=(-1{)}^m\frac{{\psi }^{(m-1)}({\Lambda }_{\alpha }/{\Delta }_{\alpha })}{{\Delta }_{\alpha }^m},$$

$$\begin{split}\rho(s,\alpha)=\frac{\Delta_{\alpha}}{\Gamma(\Lambda_{\alpha}/\Delta_{\alpha})}\exp\bigg{(}&e^{\psi^{(0)}(\Lambda_{\alpha}/\Delta_{\alpha})-\Delta_{\alpha}\,s}-\Lambda_{\alpha}s\\ &\quad+\frac{\Lambda_{\alpha}}{\Delta_{\alpha}}\psi^{(0)}(\Lambda_{\alpha}/\Delta_{\alpha})\bigg{)}.\end{split}$$

$${\tau }_2(x_i|x_f)=\frac{2}{g_0^2}{\int }_{x_i}^{x_f}dy{\int }_{-\infty }^{y}dz\tau (z|x_f)e^{-\frac{2}{g_0^2}\left[V(z)-V(y)\right]},$$

Full text

I Introduction

In this paper, we investigate deep connections between barrier crossing, the renormalization group, and the noisy saddle node bifurcation. In particular, we show that Kramers’ reaction rates can be understood as an asymptotic limit of the universal scaling near the continuous transition between high-barrier and barrier-less regimes. Applying methods from stochastic processes theory we derive an analytical expression for the universal scaling function for the mean barrier escape time near the critical point, giving the crossover between high and low barrier limits. The renormalization group provides a framework within which this result can be understood and systematically improved by perturbative calculations of corrections to scaling, some of which we give explicitly.

Barrier crossing arises in applications across physics, chemistry, and biology. In 1940, Kramers computed the barrier crossing rate for particles in both overdamped and underdamped regimes [1]. This result and others [2, 3, 4] provided the theoretical explanation for the Arrhenius equation describing chemical rate coefficients $k\sim\exp(-E_{b}/k_{B}T)$, where $E_{b}$ is the energy barrier for activation [5]. More recent efforts have established the escape rate at arbitrary damping, giving the crossover between the low- and high-damping limits [6, 7], and have accounted for the effects of state-dependent [8, 9], non-gaussian [10, 11, 12], and colored [13, 14, 12] noise, anharmonic corrections [15, 16], and fluctuating barriers [17, 14].

Most transition-state calculations assume a large barrier limit. This means the barrier escape is a rare event, with a separation of time scales between the relaxation into a quasi-equilibrium state and the escape from that state [18]. In the limit of vanishing barrier, however, there is a qualitative change in behavior. Particles instead slide down a monotonic potential, spending the most time near its inflection point. To capture the low barrier escape rate, extensions to Kramers’ theory have been developed (e.g. incorporating anharmonic corrections), but these have significant errors when the barrier and thermal energy are comparable ($E_{b}\approx k_{B}T$) [15].

Finite barrier escape problems have garnered increasing theoretical interest over the past decade, with several studies contributing further low barrier refinements of existing theories [19, 20, 21, 22, 23, 24, 20] or focusing directly on the saddle-node bifurcation where the barrier vanishes [25, 26]. Such escape processes are relevant to certain high precision measurements. For instance, force spectroscopy experiments apply a force on a single bond in a biomolecule until it breaks [27, 20]. For typical molecules, the critical force, at which the energy barrier for breaking vanishes and Kramers’ theory breaks down, is now well within the reach of atomic force microscopy and optical tweezers [27]. Another exciting application is in micro- and nano-electromechanical devices, which sensitively switch oscillation amplitude in response to an input signal by operating near the barrier-less critical point [25, 28]. Here, an analytical theory of low barrier crossing would help to distinguish between noise and signal activated switching.

We develop a critical theory for barrier crossing with a renormalization group approach that gives a complete scaling description of the noisy saddle-node bifurcation. We are inspired by previous work on the ‘intermittency’ 111In this paper ‘intermittency’ refers to the chaotic dynamical intermittency studied in Refs. [30, 31, 32] that emerges in discrete maps via a tangency bifurcation. We use this term sparingly to avoid confusion with intermittency in fluid dynamics and other areas. route to chaos [30, 31, 32], where the renormalization group coarse-grains in time, then rescales the system to fix a certain term in the potential. In chaos theory, this procedure involves iterating and rescaling a discrete map [31, 32], leading to a different fixed point for the same renormalization group equations used by Feigenbaum to study period doubling [33]. We take the continuous time limit, reducing the renormalization group to a series of elementary rescalings and yielding a simplified description applicable to barrier escape problems. Our procedure organizes what amounts to dimensional analysis, providing an elegant renormalization-group framework that unifies Kramers’ theory for Arrhenius barrier crossing with the dynamical systems theory of a noisy saddle-node bifurcation.

Why do we frame our analysis in terms of the renormalization group, if the scaling form can be justified using dimensional analysis and the analytical methods we use are drawn from more traditional stochastic analysis? On the one hand, this forms a wonderful case study, unifying and illuminating bifurcation theory, the renormalization group, and chemical reaction theory. Second, barrier crossing forms the solvable limiting case of much more complex phenomena: coupling to colored-noise heat baths, nucleation of abrupt phase transitions, and depinning transitions in disordered systems (see Section V). The scale invariance of random walks does not demand a renormalization group proof of the central limit theorem, and the one-dimensional Ising model can be solved without the machinery of flows in Hamiltonian space. But framing the problems in terms of the renormalization group provide excellent pedagogical exercises, and the natural framework for extensions to self-avoiding random walks and Ising models in higher dimensions.

Using our scaling theory as an organizing framework, we derive an analytical expression for the scaling form of the mean escape time near the saddle node bifurcation. We also compute corrections to scaling due to anharmonicity in the potential and finite initial or final position. Going beyond the mean, we develop an accurate approximation to the scaling form for the full distribution of escape times by assuming the eigenvalues of the Fokker-Plank operator are equally spaced.

As a starting point, we consider the equation of motion for an overdamped particle in a general potential $V(x)$ and driven by spatially dependent white noise,

[TABLE]

Here $f(x)=-\eta^{-1}\,dV/dx$ is the force exerted on the particle (divided by the damping coefficient $\eta$) and $g(x)$ is the spatially varying noise amplitude (with the damping absorbed). The noise $\xi(t)$ has zero mean, $\langle\xi(t)\rangle=0$ and is uncorrelated in time, $\langle\xi(t)\xi(t^{\prime})\rangle=\delta(t-t^{\prime})$. With barrier crossing phenomena in mind, we consider potentials with boundary conditions $V(x)\rightarrow\infty$ as $x\rightarrow-\infty$ and $V(x)\rightarrow-\infty$ as $x\rightarrow\infty$. The potential either has a single barrier or is monotonically decreasing (e.g. Figure 1). The quantity of interest is the mean barrier crossing time $\tau$, defined as the time particles take to reach $+\infty$ from an initial position at $-\infty$.

Besides the experimental systems discussed above, this model also serves as the natural description for a general chemical reaction, involving the transition between metastable species $A$ and $B$. These species are points in a $3N$ dimensional configuration space defined by the locations of $N$ reaction constituents. As derived by Hänggi et al. this system can be reduced to the one-dimensional model we study [18]. They coarse-grain to a one dimensional reaction coordinate, which parameterizes the minimal gradient path between the states $A$ and $B$, neglecting effects of memory friction and noise correlations, and taking the overdamped limit produces Eq. (1). The effective potential along the reaction coordinate has a barrier separating species $A$ and $B$.

II Renormalization Group and Scaling Theory

To begin, we parameterize Eq. (1) by the Taylor coefficients of $g(x)$ and $f(x)$,

[TABLE]

The renormalization group defines a flow in this space of systems described by a single reaction coordinate $x$. Near the renormalization group fixed point, the behavior is most effectively described by a single Taylor expansion at the origin. In contrast, for large barriers in Kramers’ theory, the escape time is characterized by two expansions, capturing the harmonic oscillations in the potential well and at the top of the barrier. These two equivalent schemes are shown in Fig. 1. Given the later expansion at the two extrema, the expansion at the origin can be reconstructed via a two-point Padé approximation [34].

As discussed above, the discrete renormalization group coarse-grains by iterating a map, evolving the equations forward in time. Ignoring the noise for the moment, we consider a discrete approximation to Eq. (1), $x_{n+1}=x_{n}+\delta tf(x_{n})\equiv h(x_{n})$. References [30, 31, 32] study this discrete equation using the Feigenbaum renormalization group transformation. This transformation iterates the map and rescales space, inducing a flow in function space, $T[h](x)=ah(h(x/a))$, where $a$ is the rescaling factor. The renormalization fixed point is the function $h(x)$ that obeys $T[h]=h$ for a particular rescaling $a$. We simplify this calculation by taking a continuum limit. Expanding the renormalization group to first order in $\delta t$, $T[h](x)=x+2\delta t(af(x/a))$. Thus, in the continuous time limit, the RG iteration becomes a simple rescaling of time and space. We will use this below to derive our scaling theory for barrier crossing near a saddle node bifurcation.

Within the context of the renormalization group for singular perturbations developed by Goldenfeld, Oono, and others [35, 36], our problem can be understood as having zero anomalous dimension. For our calculations in the following sections we do not need to use the traditional renormalization-group machinery; instead we assemble a variety of tools from probability theory and Markov processes, to express the mean and distributions of escape times as universal scaling functions near the transition where the barrier vanishes. The scaling theory provides a powerful and elegant structure which organizes our understanding barrier crossing.

Following the above expansion of the discrete renormalization group, we ‘coarse grain’ the system in time by scaling, $\hat{t}=t/b$. As the time-scale shrinks, the noise is amplified, $\hat{\xi}(\hat{t})=b^{1/2}\xi(t)$ (the exponent 1/2 follows from the units of the correlation function). Our goal is to understand the scaling properties near the critical point, where a qualitative change in behavior occurs. For a generic analytic potential this happens when the barrier vanishes and $V(x)=-x^{3}$ is locally a perfect cubic. Therefore, we rescale our system to fix the coefficient $\epsilon_{2}$, corresponding to the cubic term in the potential. The correct rescaling defines a new spatial coordinate $\hat{x}=bx$. After both coarse graining and rescaling, we arrive at

[TABLE]

We can then read off how the parameters flow under the renormalization group, $\hat{\epsilon}_{n}=b^{2-n}\epsilon_{n}$ and $\hat{g}_{n}=b^{3/2-n}g_{n}$. These flows and exponents exactly match those found under the discrete-time renormalization group [31, 32], indicating that the scaling of the ’intermittency route to chaos’ [29] is also non-anomalous [36]. Taking the coarse graining factor to be close to 1, $b=(1+d\ell)$, we obtain continuous flow equations,

[TABLE]

The eigenfunctions of the renormalization group in our continuum theory are the monomials $x^{n}$ and noisy monomials $\xi(t)x^{n}$. If the right hand side of Eq. (2) is an eigenfunction, it is scaled by a constant factor under the action of the renormalization group. These eigenfunctions are the much simpler continuous time limit of those for the discrete-time renormalization group [32]. In particular, the cubic potential $V(x)\propto-x^{3}$ (without noise) is the fixed point. At the fixed point, particle trajectories $x(t)\sim 1/t$ exhibit scale invariance in time as they approach the cubic inflection point at $x=0$. Perturbations away from the fixed point lead to dynamics with non-power law decay to a locally stable state or over the inflection point.

The mean barrier crossing time is a function of the potential shape and the noise correlation, encoded through the expansion coefficients $\epsilon_{n}$ and $g_{n}$. Thus, the escape time can be expressed as $\tau(\{\epsilon_{n}\},\{g_{n}\})$, where $n\geq 0$. If we coarse-grain until $g_{0}(\ell^{*})=1$, we find that the escape time has the form

[TABLE]

where $\mathcal{T}$ is a universal scaling function, with $n\geq 1$ for the second term in brackets.

While the scaling form Eq. (5) could have been written down using dimensional analysis, the renormalization group approach provides the natural structure and motivation for our approach. The parameter space flows indicate that, with a fixed quadratic force, the constant and linear force and noise terms $\{\epsilon_{0},\epsilon_{1},g_{0},g_{1}\}$ are relevant, growing under coarse graining and dominant on long time scales. Other variables are irrelevant and can be incorporated perturbatively. Of the relevant variables, the linear force coefficient $\epsilon_{1}$ can be set to zero by placing the origin at the inflection point of the potential. The spatial dependence of the noise (including the relevant linear term $g_{1}$) can also be removed by a change of coordinates $x\rightarrow\tilde{x}$ with $\tilde{x}$ defined by [37]

[TABLE]

producing a system with constant noise $\tilde{g}(\tilde{x})=g_{0}$ and force $\tilde{f}(\tilde{x})=f(\tilde{x})/g(\tilde{x})$ (hence $g_{1}$ was relevant before we removed it, because it contributes to the linear term in the expansion of $\tilde{f}$).

Systems near enough to the critical point therefore can be modeled as a cubic potential with a linear perturbation $V(x)=-x^{3}/3-\epsilon_{0}x$ and constant noise $g_{0}$. This is the ‘normal form’ used in bifurcation theory for the saddle node transition, and might have been anticipated from Taylor’s theorem. The escape time scaling form becomes,

[TABLE]

Thus, the problem asymptotically reduces to finding the universal function of a single variable $\mathcal{T}(\alpha)$, where $\alpha=\epsilon_{0}/g_{0}^{4/3}$. The limiting form of the scaling function $\mathcal{T}(\alpha)$ must give the known solutions. In the limit $\alpha\rightarrow-\infty$ the barrier is large compared to the noise, so the Kramers approximation [1] applies, so

[TABLE]

For our choice of parameters, the energy barrier is given by $E_{b}/k_{B}T=8/3|\alpha|^{3/2}$. In the opposite limit $\alpha\rightarrow\infty$, the potential is downward sloping with gradient much larger than the noise level. The passage of particles over the inflection point occurs even in the absence of noise (in contrast to the Kramers limit, which requires noise for barrier escape). Therefore, the crossing time approaches that for a deterministic particle in the cubic potential. One can easily show that the limiting scaling form is

[TABLE]

III Mean Escape Time

III.1 Analytical escape time for relevant variables

We now turn our focus to obtaining an exact analytical expression for $\mathcal{T}(\alpha)$ that is valid for all $\alpha$. To this end, we study the trajectories of particles injected at position $x_{i}$ and time $t_{i}$ into a general potential $V(x)$ with noise $g_{0}$ and compute the mean first passage time to $x_{f}$, following the standard approach [18, 38, 30]. Let $\mathcal{P}(x,t)$ be the distribution of particles over positions $x$ at time $t$, with $\mathcal{P}(x,t_{i})=\delta(x-x_{i})$. The probability that a particle has not reached $x_{f}$ at time $t$ is

[TABLE]

Note that $\mathcal{P}(0)=1$ and $\mathcal{P}(t)\rightarrow 0$ as $t\rightarrow\infty$ as long as there is noise driving the system, which guarantees particles reach $x_{f}$. The distribution of first passage times is $p(t)=-d\mathcal{P}/dt$ so that the mean first passage time is

[TABLE]

where we integrate by parts for the second equality. To derive a differential equation for $\tau(x_{i}|x_{f})$, we start from the Kolmogorov backward equation for distribution $\mathcal{P}(x,t)$ with initial condition $x_{i}$ [39],

[TABLE]

To write this equation in terms of the mean first passage time $\tau$, we multiply both sides by $t$ and integrate over $x$ and $t$. Using the relations in Eqs. (10) and (11) and the identity $d\mathcal{P}(x,t)/dt_{i}=-d\mathcal{P}(x,t)/dt$, we arrive at

[TABLE]

This gives an ordinary differential equation for the first passage time from $x_{i}$ to $x_{f}$ of particles in potential $V(x)$ and constant noise with amplitude $g_{0}$. The boundary conditions are $\tau(x_{f}|x_{f})=0$ and $\tau^{\prime}(-\infty|x_{f})=0$, which encode absorbing and reflecting boundaries respectively. Writing the solution to Eq. (13) in integral form, we arrive at the result obtained in Refs. [18, 38, 30],

[TABLE]

which satisfies the boundary conditions as long as $V^{\prime}(x)\rightarrow\infty$ as $x\rightarrow-\infty$. For large barriers, it is known that Eq. (14) reproduces Kramers escape rate formula via a saddle point approximation that expands the potential around the maximum and the minimum (as shown in Fig. 1) to second order [18].

Our renormalization group analysis allows us to restrict our focus to the relevant variables. For the cubic potential (systems on the unstable manifold of the renormalization group fixed point), the escape time can be computed analytically using Eq. (14) in the limit $x_{f}=-x_{i}\rightarrow\infty$. We find that $\tau=g_{0}^{-2/3}\mathcal{T}(\alpha)$ with the universal scaling function given by

[TABLE]

where $\text{Ai}(x)$ and $\text{Bi}(x)$ are the first and second Airy functions and $\alpha=\epsilon_{0}/g_{0}^{4/3}$ as above. This solution is shown in Figure 2, along with the Arrhenius and deterministic limits given in Eqs. (8) and (9) respectively and the mean barrier crossing times from direct simulations of the Langevin process [Eq. (1)]. The universal scaling function $\mathcal{T}(\alpha)$ reproduces the two known limits when the barrier is large or the potential is strongly downward sloping and agrees excellently with the numerical results.

Kramers’ escape rate for the cubic potential follows from Eq. (15) and the asymptotic form of the second Airy function. As $\alpha\rightarrow 0$, however, contributions from the first Airy function become important so that Kramers’ theory and extensions involving anharmonic corrections break down. The difference between Eqs. (8) and (15) is also related to the narrowing of the spectral gap of the barrier crossing Fokker-Plank operator (which has been measured numerically [40] and is discussed below in Section IV).

III.2 Corrections to scaling

III.2.1 Finite launching and absorbing positions

One limitation of our result Eq. (15) is that we assume initial and final states at infinity. In a real chemical or mechanical system, the transition of interest generally occurs between states with finite coordinates. For systems with large barriers, the escape time is exponentially large compared to the time it takes to settle into a metastable state in the well. Thus, the scaling function Eq. (15) is universal: independent of initial and final conditions. When the potential is downward sloping, the barrier crossing time is still dominated by the time spent near the inflection point and we can systematically compute corrections to scaling due to the finite initial and final positions.

Working in terms of scaling variables $\alpha=\epsilon_{/}g_{0}^{4/3}$, $\chi_{-}=x_{i}/g_{0}^{2/3}$, and $\chi_{+}=x_{f}/g_{0}^{2/3}$, we can write the barrier crossing time for the cubic potential with arbitrary initial and final conditions as $\tau=g_{0}^{-2/3}\mathcal{T}(\alpha,\chi_{-},\chi_{+})$ with the scaling function,

[TABLE]

Here $\mathcal{T}_{\pm}(\alpha,\chi_{\pm})$ are the universal corrections for finite final and initial conditions respectively. These have integral representations,

[TABLE]

and

[TABLE]

Assuming $\chi_{-}<0$ and $\chi_{+}>0$ these can be expanded in powers of $1/\chi_{\pm}$,

[TABLE]

The first term, which we can compute exactly, is simply the correction due to the deterministic trajectory in the cubic potential between $\chi_{\pm}$ and $\pm\infty$. The higher order terms describe the influence of noise on the corrections due to finite initial conditions. These terms appear to have the form $\chi_{\pm}^{-1-3j}f_{j}(\alpha/\chi_{\pm}^{2})$ for integers $j$ and some functions $f_{j}$. The expansions of $f_{1}$ and $f_{2}$ are given in parenthesis in the expression above.

Our corrections to scaling for finite launching and absorbing positions are universal if the initial conditions are sufficiently close to the cubic inflection point so that anharmonic corrections are small. Thus, for an arbitrary potential, the universal escape time near a saddle node bifurcation is the given by Eq. (15) corrected using Eq. (19). In the following section we will separately treat the anharmonic corrections to scaling. The interplay between corrections to scaling due to anharmonicity and initial conditions will be an interesting subject for future studies.

III.2.2 Anharmonic corrections

The scaling function Eq. (15) also serves as a starting point from which the theory can be systematically improved by computing anharmonic corrections to scaling. The higher order terms in the potential are irrelevant variables under the renormalization group flows and hence can be treated perturbatively. For instance, consider a quartic perturbation $\delta V(x)=-\epsilon_{3}x^{4}/4$ and let $\beta=\epsilon_{3}g_{0}^{2/3}$. In the Kramers regime, $\alpha\rightarrow-\infty$, we have that $\mathcal{T}(\alpha,\beta)\approx\mathcal{T}(\alpha)+\beta^{2}\mathcal{T}_{3}(\alpha)$ to leading order, where

[TABLE]

In the deterministic regime, $\alpha\rightarrow\infty$, we also add a quintic term as a regulator on the boundary conditions of the potential, $\delta V(x)=-\epsilon_{3}x^{4}/4-\epsilon_{4}x^{5}/5$, with $\epsilon_{4}>0$ and sufficiently large so that the potential remains monotonically decreasing. To quadratic order in $\beta$ and $\gamma=\epsilon_{4}g_{0}^{4/3}$, the universal scaling function is 222Here and in Eq. 20 we could have avoided branch cuts in our scaling functions by using the variables $\sqrt{\alpha}$, $\sqrt{\beta}$, and $\sqrt{\gamma}$. We choose our convention to avoid complex valued scaling variables and imaginary numbers in Eq. 20.

[TABLE]

The term $\pi/\sqrt{\alpha}$ is just the deterministic limit of the scaling form for the cubic potential and $\beta^{2}\mathcal{T}_{3}(\alpha)=-15\pi\sqrt{\alpha}\beta^{2}/8$ comes from the quartic perturbation to the inflection point. Other terms arise from quintic corrections or global changes in the potential. Here $\gamma$ is a dangerous irrelevant variable [42, Sections 3.6, 5.4, & 5.6], which has a pole $3\pi\beta^{2}/4\sqrt{\gamma}$ in the expansion about 0, because it is needed to keep the potential monotonic (for $\beta\neq 0)$.

IV Approximating the distribution of escape times

To completely characterize the escape times, we require their distribution, which captures the full range of outcomes we might expect from the stochastic dynamics of the Langevin equation. For high barriers (the Arrhenius limit) the escape is dominated by the decay of a single (quasi-equilibrium) mode in the bottom of the potential well, so the distribution is exponential with rate parameter given by $1/\tau$. On the other hand, for small barriers or sloping potentials many modes contribute and the mean may not be representative of the escape times in general. In this section, we develop an approximation to the distribution of escape times which is accurate for all $\alpha$ (i.e. any cubic potential). The approximate distribution is given as an analytical scaling form parameterized by two variables which are computed numerically for a given $\alpha$.

IV.1 First-passage distributions in Markov processes

To study the distribution, we borrow a result from the theory of Markov processes. Consider a birth-death process, i.e. a Markov chain with a one dimensional state space and hopping only between nearest neighbor states, that has one reflecting boundary and one absorbing boundary. This is precisely what we would obtain by discretizing space in the Langevin equation Eq. (1), which has a one dimensional state space with a reflecting boundary for large negative $x$ values and an absorbing boundary for large positive $x$ values. For Markov systems initialized near the reflecting boundary, the distribution of times to reach the absorbing state can be characterized completely in terms of the eigenvalues of the transition matrix. In particular the distribution, $p(t)=-d\mathcal{P}(t)/dt$ (using the notation from the preceding section) can be shown to be a convolution of exponentials [43],

[TABLE]

where $\lambda_{n}$ are the negative eigenvalues of the Markov transition matrix, $\mathcal{E}(\lambda_{n})$ are exponential distributions, and $*$ denotes a convolution. As we will see below, it is useful consider the Fourier transform $\tilde{p}(\omega)$ of the distribution $p(t)$ and the cumulant generating function

[TABLE]

since the cumulants $\kappa_{m}$ are easily expressed in terms of the eigenvalues of the Markov matrix,

[TABLE]

Note that the first cumulant $\kappa_{1}$ is just the mean barrier crossing time $\tau$, which was the focus of the previous section.

IV.2 Spectra of the Fokker-Planck operator

When the state space becomes continuous, the process is generated by the Fokker-Planck operator rather than a finite transition matrix. Thus we want to understand the spectra of the right hand side of the Fokker-Planck equation,

[TABLE]

For numerical evaluation of the eigenvalues is convenient to change variables so that the differential operator is Hermitian. To do this we rescale $\mathcal{P}(x,t)$ by the square root of the Boltzmann factor, defining $\sigma(x,t)=\mathcal{P}(x,t)/\exp(-V(x)/g_{0}^{2})$. Then $\sigma(x,t)$ satisfies

[TABLE]

which has the form of a Schrödinger equation with imaginary time. Importantly, even in the limit of continuous space, the spectrum of the Fokker-Plank operator is discrete. This is easy to see from Eq. (26), because the ‘effective quantum potential’ $(V^{\prime}(x))^{2}/2g_{0}^{2}-V^{\prime\prime}(x)/2$ is bounded from below for the cubic potential as well as any polynomial potential that diverges super-linearly as $x\rightarrow\pm\infty$.

Within the scaling theory developed in Section II, if we specialize to the cubic potential (i.e. relevant variables of the renormalization group) the eigenvalues have scaling forms

[TABLE]

written in terms of the RG-invariant quantity $\alpha=\epsilon/g_{0}^{4/3}$. This also implies a scaling form for the cumulants $\kappa_{m}=g_{0}^{-2m/3}\mathcal{K}_{m}(\alpha)$. The scaling forms $\Lambda_{n}(\alpha)$ are plotted in Fig. (3) for the first several eigenvalues. These curves describing the $\alpha$-dependence of the Fokker-Planck spectra are universal for systems near a saddle node bifurcation (where $\alpha$ is the only relevant variable). For large positive $\alpha$, the eigenvalues are approximately evenly spaced, while for large negative $\alpha$, a single slowly decaying mode with eigenvalue very close to 0 dominates the behavior.

IV.3 Evenly spaced eigenvalue approximation

To develop an approximation to the distribution of barrier escape times for the cubic potential, we assume the eigenvalues are equally spaced, $\Lambda_{n}(\alpha)=\Lambda_{\alpha}+n\Delta_{\alpha}$, where $\Lambda_{\alpha}\equiv\Lambda_{0}(\alpha)$ is the decay rate of the slowest decaying eigenmode. This approximation becomes exact in the limit $\alpha\rightarrow\infty$, for which the effective potential in Eq. (26) becomes harmonic. Though we can see in Fig. (3) that this approximation is clearly not correct for $\alpha\rightarrow-\infty$, the behavior in this limit is dominated by the vanishing leading eigenvalue $\Lambda_{\alpha}$. For instance, the cumulants in Eq. (24) are insensitive to $\Delta_{\alpha}$ when $\Delta_{\alpha}\gg\Lambda_{\alpha}$ and $\Lambda_{\alpha}\rightarrow 0$. Thus, this approximation will produce a family of distributions which correctly captures the behavior for both large negative and positive $\alpha$. As we will see below, the approximation is also quite accurate over the full range of $\alpha$.

One caveat of the equal spacing approximation is that the mean escape time $\mathcal{T}=\mathcal{K}_{0}=\sum_{n=1}^{\infty}\Lambda_{n}^{-1}$ always diverges. Therefore, we must fix the mean to the value derived in Section III.1. On the other hand, the higher order cumulants, which determine the shape of the distribution, always converge.

To proceed, we sum Eq. (24) using our eigenvalue ansatz. For $m\geq 2$ the result is,

[TABLE]

where $\psi^{(m)}(x)$ is the polygamma function of order $m$. Neglecting the mean, we can sum the series Eq. (23), then exponentiate and inverse Fourier transform to obtain the distribution $p(t)$. Writing the distribution in terms of its universal scaling function, $p(t,\epsilon,g_{0})=g_{0}^{2/3}\rho(s,\alpha)$ with $s=(t-\tau)/g_{0}^{2/3}$ and $\alpha=\epsilon/g_{0}^{4/3}$, the final result is

[TABLE]

To use this family of distributions, we compute the leading eigenvalue $\Lambda_{\alpha}$ using a shooting method on the right hand side of Eq. (26). The effective $\Delta_{\alpha}$ is chosen so that the second cumulant $\mathcal{K}_{2}$ (i.e. the variance) is exact. To do this, we follow the same approach used in Section III.1 to write an integral expression for the variance. The result is $\mathcal{K}_{2}=\tau_{2}-\tau^{2}$ with

[TABLE]

where $\tau(z|x_{f})$ is the mean given in Eq. (14). After evaluating these integrals numerically to obtain the variance, we solve $\mathcal{K}_{2}=\psi^{(1)}(\Lambda_{\alpha})/\Delta_{\alpha}^{2}$ to fix the eigenvalue spacing $\Delta_{\alpha}$.

In Figure 4 we compare the approximate distribution Eq. (29) to those obtained by direct numerical simulation of the Fokker-Plank equation. Our approximation captures the shape of the distributions remarkably well for the full range of $\alpha$.

We can also show analytically that this distribution reproduces the correct large $\alpha$ limits. For large negative $\alpha$, taking the $\Delta_{\alpha}\rightarrow\infty$ limit of Eq. (29), gives the expected exponential distribution, $\rho(s,\alpha)\to\Lambda_{\alpha}e^{-s\Lambda_{\alpha}-1}$ (shifted to have zero mean). For large positive we can estimate $\Lambda_{\alpha}$ and $\Delta_{\alpha}$ using the harmonic approximation to the potential in Eq. (26). We find $\Lambda_{\alpha}\sim\alpha^{2}$ and $\Delta_{\alpha}\sim\sqrt{\alpha}$ so that the cumulants scale like $\mathcal{K}_{m}(\alpha)\sim\alpha^{3(1-m)/2}$ for $\alpha\to\infty$. In particular, the ratio $\mathcal{K}_{m}/\mathcal{K}_{2}^{m/2}\to 0$ for $m>2$, i.e. the higher order cumulants are small compared to the variance and the distribution approaches a Gaussian for large $\alpha$. This prediction is born out in direct simulation of the Fokker-Plank equation (see Figure 4d).

The distribution Eq. (29) combined with our analytical understanding of the mean escape time provides a complete description of the barrier crossing process. Our procedure allows for accurate approximations to the distribution by evaluating just the two universal scaling forms $\Lambda_{\alpha}$ and $\Delta_{\alpha}$. Knowledge of these quantities allows for evaluation of the distribution of barrier crossing times for any system near a saddle-node bifurcation.

V Discussion

We expect our results will be directly applicable to barrier crossing processes in which thermal fluctuations are comparable to the energy barrier including the aforementioned experimental systems, narrow escape problems in cellular biology [44], and downhill protein folding scenarios [45, 46]. Our approximation to the distribution of barrier escape times, combined with our analytical results for the mean, provides an accurate and complete characterization of the barrier crossing process for systems near a saddle node bifurcation.

A more thorough analysis of incorporating perturbative corrections from irrelevant variables into Eq. (14) would be both theoretically interesting and useful in applications. The interplay between anharmonic corrections and finite initial conditions is also important. Computing these corrections will extend the applicability of our theory to systems in which the boundary conditions correspond to positions with non-negligible anharmonicity. It will be interesting to test the accuracy of our evenly spaced eigenvalue approximation to the barrier escape time distribution when irrelevant variables are incorporated. We conjecture that, at least perturbatively, Eq. (29) will still accurately parameterize the distributions if $\Lambda_{\alpha}$ and $\Delta_{\alpha}$ are corrected to account for the irrelevant variables.

Our analysis directly translates to higher order cuspoid catastrophes [47, Section 36.2], which form their own universality classes with different exponents (in fact, these have already been analyzed for the discrete iterated map [31, 32]). For these bifurcations, the fixed point potential will be a higher order monomial and our analysis can be used to identify the relevant variables and develop a scaling theory for quantities like the barrier crossing time mean and distribution.

More generally, it would be useful to study the applicability of our renormalization group and scaling analysis to systems with colored noise, multiple dimensions, or in other damping regimes. The effects of colored noise are encoded in the correlation function $\langle\xi(t)\xi(t^{\prime})\rangle=G(x,t-t^{\prime})$. The renormalization group transformation can be adapted to act on the Fourier transform of this quantity $\tilde{G}(x,\omega)$, giving flows of the colored noise under coarse-graining (for example, barrier crossing between two symmetric wells – a noisy pitchfork bifurcation – coupled to an Ohmic heat bath leads in the quantum limit to a critical point in the same universality class as the Kondo problem [48]). We expect short-range correlations will be irrelevant under coarse-graining, while those with power-law decay will give rise to new anomalous scaling. For some reactions, an underdamped model or multi-dimensional reaction coordinate may be required for an accurate description. Renormalization group scaling will provide a natural organizing framework for these studies.

Nucleation of abrupt phase transitions (e.g., raindrop formation) is also described by Arrhenius rates. Here the noiseless bifurcation underlying critical droplet theory is the spinodal line. This line – a discredited mean-field boundary between nucleation and spontaneous phase separation – could play the role of our renormalization-group fixed point in a future generalization of this work to higher spatial dimensions.

Finally, the saddle node bifurcation is the simplest example of a depinning transition. One anticipates studying how adding noise would affect depinning of earthquakes, vortices in superconductors, plastic flow in crystals, raindrops on windshields, coffee soaking into napkins, and other depinning phenomena. Each of these has anomalous exponents and RG treatments [49] even without added noise. Our work provides a stepping stone toward understanding the universal scaling near noisy depinning transitions in these more sophisticated systems as well.