First Passage Times: the Kramers problem

This tutorial illustrates the use of the stochrare package for first-passage time computations in the context of the Kramers problem [1] (diffusion in a double well potential). Here, we simply compute numerically the rate of transition between the two attractors.

References

• [1] Kramers, Physica, 7, 284-304 (1940)
• [2] Gardiner, Handbook of Stochastic Methods, Springer, chap. 9 and §5.5.
• [3] Risken, The Fokker-Planck equation, Springer, §5.10
• [4] Caroli, Caroli and Roulet, J. Stat. Phys., 21, 415-437 (1979)
• [5] Caroli, Caroli and Roulet, J. Stat. Phys., 26, 83-111 (1981)
• [6] Hanggi, Talkner, Borkovec, Rev. Mod. Phys., 62, 251-342 (1990)

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import stochrare as sr

Let us define the DoubleWell class, corresponding to a simple bistable system, by subclassing the generic 1D diffusion class ConstantDiffusionProcess1D:

class DoubleWell(sr.dynamics.diffusion1d.ConstantDiffusionProcess1D):
    """
    Double well potential model.
    """
    default_dt = 0.01

    def __init__(self, Damp, **kwargs):
        sr.dynamics.diffusion1d.ConstantDiffusionProcess1D.__init__(self, lambda x, t: -x*(x**2-1), Damp, **kwargs)

    def potential(self, X, t):
        """
        Return the value of the potential at the input points.
        """
        Y = X**2
        return Y*(Y-2.0)/4.0

The dynamics is given by the stochastic differential equation:

\[ \begin{align}\begin{aligned}dX_t = -V'(X_t)dt+\sqrt{2D}dW_t,\\with :math:`V(x)=x^4/4-x^2/2`.\end{aligned}\end{align} \]

Qualitative understanding

Deterministic dynamics

Let us first look at the dynamics for the deterministic system: we plot below the phase portrait of the system, superimposed upon the potential.

ax = plt.axes(xlabel='t', ylabel='x(t)')
X = np.linspace(-1.5, 1.5, num=100)
ax.contourf(np.linspace(0, 10), X, np.tile(DoubleWell(0).potential(X, 0), (50, 1)).T, 50, cmap='copper')
for x0 in X:
    ax.plot(*DoubleWell(0).trajectory(x0, 0, T=10), color='white')

\(x=-1\) and \(x=1\) are stable fixed points with basins of attraction \(]-\infty,0[\) and \(]0,+\infty[\), respectively, separated by an unstable fixed point at \(x=0\). Typical relaxation time is of order one but diverges as \(x_0\) goes to 0.

Noise-induced transitions between the two attractors

When the diffusion coefficient does not vanish, one should expect transitions between the two attractors. Let us illustrate this below with a moderate value of the noise.

DoubleWell.trajectoryplot(DoubleWell(0.1, deterministic=True).trajectory(-1, 0, T=1000, dt=0.01));

The stationary distribution of the system should therefore be bimodal. It is simply given by \(e^{-V(x)/D}\), up to a normalization factor. Let us compute this stationary distribution numerically by solving the Fokker-Planck equation with absorbing boundary conditions at a sufficient distance.

_, ax = DoubleWell(0.1).pdfplot(0.0, 10.0, 100.0, dt=0.0005, npts=600, bounds=(-3.0, 3.0),
                                P0=sr.fokkerplanck.FokkerPlanck1D.gaussian1d(-1, 0.1, np.linspace(-3, 3, 600)));
x = np.linspace(-3, 3)
V = DoubleWell(0.1).potential(x, 0)
ax.plot(x, np.exp(-V/0.1)/np.trapz(np.exp(-V/0.1), x=x), color='black', ls='dotted', lw=2);

/Users/corentin/codes/stochrare/stochrare/fokkerplanck.py:151: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if P0 == 'gauss':
/Users/corentin/codes/stochrare/stochrare/fokkerplanck.py:153: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if P0 == 'dirac':
/Users/corentin/codes/stochrare/stochrare/fokkerplanck.py:157: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if P0 == 'uniform':

Particles initially concentrated in the left well tunnel through the potential barrier, and ultimately the relative probability of the two attractors is the same (because the quasi-potential is symmetric).

To study the statistical properties of the transitions between the two attractors, we introduce the first-passage time:

\[ \begin{align}\begin{aligned}\tau_M = \inf \{ t>0, x_t > M | x_0=-1 \}.\\It is a random variable entirely determined by the stochastic process :math:`X`. In ``stochrare``, it can be represented by the ``FirstPassageProcess`` class:\end{aligned}\end{align} \]

tau = sr.firstpassage.FirstPassageProcess(DoubleWell(0.1))

Methods of the FirstPassageProcess class allow for sampling the random variable, estimating its mean, moments and PDF, using direct simulation or using the Fokker-Planck equation.

tau_samples = tau.escapetime_sample(-1, 0, 0, ntraj=10000)

tau.escapetime_pdfplot(tau.escapetime_pdf(tau_samples))

When the noise is small enough, it is expected that the transition times are Poisson distributed. Their statistics is entirely determined by the transition rate (the parameter of the Poisson distribution), which is the inverse of the average time between two transitions.

ax = plt.axes(xlabel=r'$\tau_M$', ylabel=r'$p(\tau_M)$', yscale='log', ylim=(1e-6, 0.1))
ax.scatter(*tau.escapetime_pdf(tau_samples), label='MC');
t = np.linspace(*ax.get_xlim())
ax.plot(t, np.exp(-t/np.mean(tau_samples))/np.mean(tau_samples), color='C1', label='Poisson');
ax.legend();

Computing the transition rate

Let us now compute the transition rate, or equivalently, the mean first-passage time. For 1D homogeneous processes, a theoretical result can be obtained analytically:

\[ \begin{align}\begin{aligned}\mathbb{E}[\tau_M] = \frac{1}{D} \int_{-1}^{M} dx e^{V(x)/D} \int_{-\infty}^x e^{-V(y/D}dy.\\In the small noise limit :math:`D \to 0`, a saddle-point approximation yields the Eyring-Kramers formula:\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned}\mathbb{E}[\tau_M] \approx \frac{2\pi}{\sqrt{|V''(0)V''(-1)|}}e^{\Delta V/D},\\for :math:`M> 0` and :math:`\Delta V=V(0)-V(-1)=1/4` here.\end{aligned}\end{align} \]

mfpt_mc = np.array([sr.firstpassage.FirstPassageProcess(DoubleWell(D)).escapetime_avg(-1, 0, 0.5, ntraj=1000)
                    for D in 1./np.arange(1, 11)])

ax = plt.axes(xlabel=r'$1/D$', ylabel=r'$\mathbb{E}[\tau_M]$', yscale='log', ylim=(1, 100))
ax.scatter(np.arange(1, 11), mfpt_mc, label='MC');
ax.plot(np.linspace(1, 10), np.sqrt(2)*np.pi*np.exp(0.25*np.linspace(1, 10)), label='Theory', color='C1')
ax.legend();

The code also allows for computing the mean first passage time using the Fokker-Planck equation or its adjoint. We will soon update this notebook to show this method.

One could also show how the prefactor of the first-passage time depends on the threshold \(M\) (above we chose \(M=0.5\)), for a fixed noise amplitude \(D\). For \(D\) small enough, there is a sharp transition around \(M=0\), and the first-passage time depends very little on \(M\) away from this boundary layer (and it is given by the above Eyring-Kramers formula). In general, the expression of the first-passage time is always given by the integral formula before the saddle-point approximation.

Finally, these transitions can also be characterized using the instanton formalism. We shall also illustrate this in a future version of the notebook.