# Weak-noise limit of a piecewise-smooth stochastic differential equation

Yaming Chen<EMAIL_ADDRESS>School of Mathematical Sciences, Queen
Mary University of London, London E1 4NS, United Kingdom  Adrian Baule
<EMAIL_ADDRESS>School of Mathematical Sciences, Queen Mary University of
London, London E1 4NS, United Kingdom  Hugo Touchette<EMAIL_ADDRESS>National Institute for Theoretical Physics, Stellenbosch 7600, South Africa
Institute of Theoretical Physics, University of Stellenbosch, Stellenbosch
7600, South Africa  Wolfram Just<EMAIL_ADDRESS>School of Mathematical
Sciences, Queen Mary University of London, London E1 4NS, United Kingdom

(October 11, 2013)

###### Abstract

We investigate the validity and accuracy of weak-noise (saddle-point or
instanton) approximations for piecewise-smooth stochastic differential
equations (SDEs), taking as an illustrative example a piecewise-constant SDE,
which serves as a simple model of Brownian motion with solid friction. For
this model, we show that the weak-noise approximation of the path integral
correctly reproduces the known propagator of the SDE at lowest order in the
noise power, as well as the main features of the exact propagator with higher-
order corrections, provided the singularity of the path integral associated
with the nonsmooth SDE is treated with some heuristics. We also show that, as
in the case of smooth SDEs, the deterministic paths of the noiseless system
correctly describe the behavior of the nonsmooth SDE in the low-noise limit.
Finally, we consider a smooth regularization of the piecewise-constant SDE and
study to what extent this regularization can rectify some of the problems
encountered when dealing with discontinuous drifts and singularities in SDEs.

###### pacs:

02.50.–r, 05.40.–a, 46.55.+d, 46.65.+g

##  I Introduction

The study of piecewise-smooth differential equations is a relatively recent
topic in field dynamical systems Filippov (1988); di Bernardo et al. (2008);
Makarenkov and Lamb (2012); Biemond et al. (2012); Colombo et al. (2012),
despite their wide use in applied sciences and engineering. The dynamics
generated by these equations displays many unexpected phenomena, including
stick-slip transitions associated, for example, with solid friction forces
Bowden and Tabor (1950); Elmer (1997); Persson (1998); Berger (2002) and
bifurcations that do not appear in the standard classification of catastrophes
of smooth dynamical systems Makarenkov and Lamb (2012); Biemond et al. (2012);
Colombo et al. (2012). They also show, in the case of systems with
discontinuous derivatives or forces (so-called Filippov systems Filippov
(1988)), multivalued solutions for a given initial condition, leading to a
loss of determinism Jeffrey (2011).

Stochastic versions of piecewise-smooth dynamical systems perturbed by noises
are also used as models of physical and biological systems. Stochastic
differential equations (SDEs) with piecewise-smooth drifts are commonly used,
for example, in stochastic or Brownian ratchets, which serve as models of
diffusion and transport in a variety of biological motors Reimann (2002).
Another important class of problem concerns diffusion of solid objects on
solid surfaces, which can be modeled phenomenologically using piecewise-smooth
SDEs with solid friction forces Caughey and Dienes (1961); Atkinson and
Caughey (1968); de Gennes (2005); Hayakawa (2005). The dynamics in this case
shows stick-slip transitions, as in the noiseless case, but also new features
due to the noise such as directed motion in the absence of a mean force bias
Baule and Sollich (2012) and noise-dependent decay of correlation functions
Touchette et al. (2010). Some of these features have been investigated
experimentally in Chaudhury and Mettu (2008); Goohpattader et al. (2009);
Goohpattader and Chaudhury (2010); Gnoli et al. (2013).

The theory of piecewise-smooth SDEs is only in its infancy compared to its
noiseless counterpart. Exact solutions for the transition probability
distribution or propagator of a few simple piecewise-constant or piecewise-
linear SDEs are known (see, e.g., Caughey and Dienes (1961); Touchette et al.
(2010); Simpson and Kuske (2012); Karatzas and Shreve (1984); Zhang (1990);
Benes et al. (1980); Browne and Whitt (1995)). Some studies have also looked
at the large deviations of SDEs with discontinuous drift (so-called SDEs with
discontinuous statistics) Korostelev and Leonov (1993); Chiang and Sheu
(2000); Boué et al. (2000); Gradinaru et al. (2001). However, from these
disconnected studies it is not clear how nonsmooth SDEs can be studied with
techniques developed and used for smooth systems. Here we focus on two such
techniques, namely, the path integral representation of propagators and the
weak-noise (also called saddle-node or instanton) approximation of these
integrals Chaichian and Demichev (2001); Graham (1989); Touchette (2009). The
weak-noise limit is particularly interesting from a physical point of view
because a piecewise-smooth system does not necessarily behave continuously
with the magnitude of a force or noise and may therefore behave in a non-
trivial way in the limit of vanishing noise. This raises the question of the
validity of the weak-noise limit for nonsmooth systems.

In this paper, we address this question in a practical way by studying a
simple piecewise-constant SDE, defined by

$\dot{v}(t)=-\mu\sigma(v(t))+\sqrt{D}\,\xi(t),$ (1)

where $\sigma(v)$ denotes the sign of the state variable $v(t)$, $\mu$ is a
positive constant, and $D$ is the strength of the Gaussian white noise
$\xi(t)$ characterized by

$\langle\xi(t)\rangle=0,\quad\langle\xi(t)\xi(t^{\prime})\rangle=2\delta(t-t^{\prime}).$
(2)

The notation $\langle\cdots\rangle$ stands for the average over the noise.
Physically, the SDE (1) is the simplest model of stick-slip motion where
$v(t)$ represents the velocity of a solid object of unit mass sliding over a
surface with solid (also called dry or Coulomb) friction coefficient $\mu$ per
unit mass Caughey and Dienes (1961); Atkinson and Caughey (1968); de Gennes
(2005); Hayakawa (2005). Mathematically, it is also the simplest piecewise-
smooth SDE whose time-dependent propagator $p(v,t|v_{0},0)$, representing the
probability density that $v(t)=v$ given the initial condition $v(0)=v_{0}$, is
known exactly Touchette et al. (2010). It is thus a good starting point for
benchmarking results about nonsmooth SDEs.

In the present case, we have

$p(v_{t},t|v_{0},0)=\frac{\mu}{D}\,\hat{p}\left(\frac{\mu}{D}v_{t},\frac{\mu^{2}}{D}t\bigg{|}\frac{\mu}{D}v_{0},0\right),$
(3)

where

$\hat{p}(x,\tau|x^{\prime},0)=\frac{e^{-\tau/4}}{2\sqrt{\pi\tau}}\,e^{-(|x|-|x^{\prime}|)/2}\,e^{-(x-x^{\prime})^{2}/(4\tau)}+\frac{e^{-|x|}}{4}\bigg{[}1+\mathrm{erf}\left(\frac{\tau-(|x|+|x^{\prime}|)}{2\sqrt{\tau}}\right)\bigg{]}$
(4)

is the propagator in nondimensional units and $\mathrm{erf}(x)$ is the error
function. This propagator is obtained by explicitly solving the corresponding
time-dependent Fokker-Planck equation, which in nondimensional units has the
form

$\frac{\partial}{\partial\tau}\hat{p}(x,\tau|x^{\prime},0)=\frac{\partial}{\partial
x}[\sigma(x)\hat{p}(x,\tau|x^{\prime},0)]+\frac{\partial^{2}}{\partial
x^{2}}\hat{p}(x,\tau|x^{\prime},0).$ (5)

Taking the limit $t\rightarrow\infty$ of $p(v,t|v_{0},0)$, we obtain the
stationary probability density function of Eq. (1),

$p(v)=\frac{\mu}{2D}\,e^{-\mu|v|/D},$ (6)

which solves the time-independent Fokker-Planck equation Hayakawa (2005). Note
that $p(v_{t},t|v_{0},0)$ and, by extension, $p(v)$ are symmetric under the
change $v_{0}\rightarrow-v_{0}$ and $v(t)\rightarrow-v(t)$ due to the
symmetric force $\sigma(v)$. Thus we can confine the analysis to the case
$v_{0}>0$ without loss of generality.

In the next section, we compare these exact results, which are valid for any
noise power $D$, with various low-noise approximations of the path integral
representation of $p(v,t|v_{0},0)$ in order to test the validity and accuracy
of these approximations and to discuss subtle singularities arising in the
path integral when dealing with discontinuous drifts. Related path integrals
were studied in Baule et al. (2010, 2011) for a non-exactly solvable model of
solid friction. For the model that we consider, we will see that the weak-
noise approximation gives the correct propagator at the lowest order in $D$,
as well as its main features with higher-order corrections, provided some
heuristics are used to treat the singularities of the path integral. The
higher-order corrections are able, in particular, to reproduce the short- and
long-time behavior of the propagator, as well as its tail behavior.

To treat the singularity of the system in a more explicit and systematic way,
we then consider in Sec. III a regularized version of our model, given by

$\dot{v}(t)=-\mu\tanh\left(\frac{v(t)}{\varepsilon}\right)+\sqrt{D}\,\xi(t),$
(7)

which recovers the original piecewise-smooth model of Eq. (1) in the limit
$\varepsilon\rightarrow 0$. Although we do not have an explicit expression for
the propagator of this smooth (but nonlinear) model, we show with Langevin
simulations that the weak-noise approximation of the path integral, which is
now well-defined and shows no singularity, reproduces the main features of the
propagator at different orders of approximation, with roughly the same
accuracy as for the singular model. The regularized model also allows us to
obtain analytical results about the instanton or optimal path of the system,
which is the most probable path singled out by the weak-noise approximation,
and the so-called action functional or quasipotential, obtained by
approximating the propagator at the lowest order in $D$ with the optimal path
Graham (1989); Touchette (2009); Freidlin and Wentzell (1984). From these
results, presented in Sec. IV, we are able to study the evolution of the
propagator’s tails in the long-time limit and to understand the behavior of
the SDE in the noiseless limit. Conclusions drawn from these results are given
in Sec. V.

##  II Piecewise-smooth model

In this section, we compare the exact propagator of the piecewise-constant SDE
of Eq. (1) with various approximations of the path integral representation of
this propagator so as to discuss the validity of these approximations for a
nonsmooth SDE. The path integral has the form

$p(v_{t},t|v_{0},0)=\int_{(v_{0},0)}^{(v_{t},t)}\mathcal{D}[v]\,J[v]\,e^{-S^{(0)}[v]/(4D)}$
(8)

and involves two terms: the action functional

$S^{(0)}[v]=\int_{0}^{t}[\dot{v}(s)+\mu\sigma(v(s))]^{2}\,\mathrm{d}s,$ (9)

which is a measure of the probability of a path $\\{v(s)\\}_{s=0}^{t}$ in
velocity space, and the Jacobian functional $J[v]$, which is the Jacobian of
the transformation of the Gaussian white noise $\xi(t)$ to the $v(t)$ process
(see Appendix A).

For nonsmooth SDEs, two problems arise with the path integral (8). The first
is that, since the noiseless system $\dot{v}=-\mu\sigma(v)$ admits in general
piecewise-linear trajectories that are continuous but nondifferentiable at
points where $v(t)$ vanishes, the minimization of the action must also be
carried over these trajectories, which means that care must be taken with the
$\dot{v}$ term in $S^{(0)}[v]$. The second problem is that the Jacobian, as
given in Eq. (70) of Appendix A, is singular. For our model, we formally have

$J[v]=\exp\left(\mu\int_{0}^{t}\delta(v(s))\,\mathrm{d}s\right)\,,$ (10)

where $\delta(v)$ is the Dirac delta function. Below we show how to treat this
singular contribution and how its inclusion or noninclusion in the saddle-
point approximation of the path integral determines different orders of
approximation of the propagator as $D\rightarrow 0$.

###  II.1 Zeroth-order saddle-point approximation

The lowest order approximation of the propagator $p(v,t|v_{0},0)$ is obtained
in the noiseless limit $D\rightarrow 0$ by finding the path
$\\{v^{(0)}(s)\\}_{s=0}^{t}$ that minimizes the action $S^{(0)}[v]$ so as to
write

$p^{(0)}(v,t|v_{0},0)=Ne^{-S^{(0)}[v^{(0)}]/(4D)},$ (11)

where $N$ is a normalization constant. We refer to this approximation as the
zeroth-order saddle-point approximation [SPA(0)]. The rationale for this
approximation is that, in the limit $D\rightarrow 0$, the path integral (8) is
dominated by the probability of the most probable or optimal path
$\\{v^{(0)}(s)\\}_{s=0}^{t}$ having minimal action. The Jacobian $J[v]$ can be
neglected at this level of approximation, since it does not depend on $D$.

The Euler-Lagrange equation associated with the minimization of $S^{(0)}[v]$
in the region $v>0$ and $v<0$ is simply

$\ddot{v}^{(0)}=0$ (12)

and leads to straight paths

$v^{(0)}_{d}(s)=(v_{0}-v_{t})s/t+v_{0},$ (13)

which we call direct paths. As mentioned above, in addition to these paths, we
must consider paths that follow the attractor $v=0$, since these paths appear
in the noiseless system. As a result, the minimization of $S^{(0)}[v]$ must be
carried out over all continuous and piecewise-linear paths consisting of
direct paths and paths following the $v=0$ axis. This situation differs from
smooth SDEs, for which the minimization is generally over all continuously
differentiable paths, and leads us to define two important heuristic
principles for dealing with nonsmooth SDEs: (i) The action of a path must be
evaluated as the sum of the actions of all its linear (or in general smooth)
parts without regard to its joining (nonsmooth) points and (ii) any part of a
path on the $v=0$ axis (or, in general, on an attractor of the noiseless
system) must have a zero action, in analogy with smooth systems.

Figure 1: (Color online) Different paths considered for minimizing the action,
Eq. (9), for (a) $v_{t}>0$ and (b) $v_{t}<0$. Here $v_{d}^{(0)}$ denotes the
direct path (13), $v_{ind}^{(0)}$ the indirect path (15), and $v_{int}^{(0)}$
the intermediate path (15) with $t_{1}=t_{2}=t_{0}$.

In the present model, two types of optimal paths arise from the action
minimization. The first consists of direct paths $v_{d}^{(0)}(s)$, found
above, which directly link the positive initial velocity $v_{0}$ to a final
velocity $v_{t}$ [see Fig. 1(a)] and whose action is

$S^{(0)}[v^{(0)}_{d}]=\left\\{\begin{array}[]{lll}(v_{t}-v_{0}+\mu
t)^{2}/t-4\mu v_{t}&&\mbox{for }v_{t}<0\\\ (v_{t}-v_{0}+\mu
t)^{2}/t&&\mbox{for }v_{t}>0.\end{array}\right.$ (14)

The second type of path is the piecewise linear path mentioned above,
consisting of two straight lines in the region $v>0$ or $v<0$ connected by a
straight line at $v=0$. The equation of these indirect paths is

$v_{ind}^{(0)}(s,t_{1},t_{2})=\left\\{\begin{array}[]{lll}(t_{1}-s)v_{0}/t_{1}&&\mbox{for
}s<t_{1}\\\ 0&&\mbox{for }t_{1}<s<t_{2}\\\ (s-t_{2})v_{t}/(t-t_{2})&&\mbox{for
}s>t_{2},\\\ \end{array}\right.$ (15)

where $t_{1}<t_{2}$ are arbitrary times at which a path reaches $v=0$ (see
Fig. 1). Following the two principles above, we evaluate the action of these
paths in a piecewise way with $S^{(0)}[0]=0$ and minimize it for $t_{1}<t_{2}$
to obtain

$S^{(0)}[v_{ind}^{(0)}]=4\mu|v_{t}|$ (16)

for $t_{1}=v_{0}/\mu$ and $t_{2}=t-|v_{t}|/\mu$. Since we require
$t_{1}<t_{2}$, the lower bound is reached if $|v_{t}|<\mu t-v_{0}$. Note that
paths arising in the limit case where $t_{1}=t_{2}\equiv t_{0}$ with
$v_{0}v_{t}<0$ correspond to direct paths, whereas those that just “bounce” on
the $v=0$ axis, i.e., $t_{1}=t_{2}\equiv t_{0}$ but $v_{0}v_{t}>0$, are called
intermediate paths [see Fig. 1(a)] and have an action equal to

$\displaystyle S^{(0)}[v_{int}^{(0)}]=(|v_{t}|+v_{0}-\mu t)^{2}/t+4\mu|v_{t}|$
(17)

for $t_{0}=v_{0}t/(v_{0}+|v_{t}|)$.

It can be easily checked that any piecewise-linear paths other than those
considered above have a greater action, and so cannot be optimal. Therefore,
combining Eqs. (14), (16) and (17) we can write the equation of the optimal
path as

$v^{(0)}(s)=\left\\{\begin{array}[]{lll}v_{d}^{(0)}(s)&&\mbox{for
}{\color[rgb]{0,0,0}t<v_{0}/\mu}\\\
v_{ind}^{(0)}(s,v_{0}/\mu,t-|v_{t}|/\mu)&&\mbox{for
}{\color[rgb]{0,0,0}t>v_{0}/\mu}\mbox{ and }v_{t}\in[v^{-}(t),v^{+}(t)]\\\
v_{d}^{(0)}(s)&&\mbox{for }{\color[rgb]{0,0,0}t>v_{0}/\mu}\mbox{ and
}v_{t}\not\in[v^{-}(t),v^{+}(t)],\end{array}\right.$ (18)

where the limits of the velocity interval are defined by

$v^{-}(t)=v_{0}-\mu t,\qquad v^{+}(t)=\big{(}\sqrt{v_{0}}-\sqrt{\mu
t}\big{)}^{2}.$ (19)

This shows that, if the endpoint $(v_{t},t)$ lies in the area bounded by
$v^{-}(s)$ and $v^{+}(s)$ with $s>v_{0}/\mu$, as shown in Fig. 2, then the
optimal path is an indirect path, otherwise it is a direct path. Intermediate
paths do not enter in this result, but will be useful when we treat the
Jacobian.

Figure 2: (Color online) Regions in the $v_{t}$-$t$ plane for which the
optimal path of the action at leading order (9) is given by an indirect
(direct) path [shaded (white)] [see Eqs. (19)]. Here $v_{c}=\mu/v_{0}$ denotes
the tip of the region.

The SPA(0) of the propagator is obtained from this result by substituting the
corresponding action in (11). For $t<v_{0}/\mu$ we find

$p^{(0)}(v_{t},t|v_{0},0)=N_{1}\left\\{\begin{array}[]{lll}\exp[-(v_{t}-v_{0}+\mu
t)^{2}/(4Dt)-\mu v_{t}/D]&&\mbox{for }v_{t}<0\\\ \exp[-(v_{t}-v_{0}+\mu
t)^{2}/(4Dt)]&&\mbox{for }v_{t}>0,\end{array}\right.$ (20)

whereas for $t>v_{0}/\mu$ we find

$p^{(0)}(v_{t},t|v_{0},0)=N_{2}\left\\{\begin{array}[]{lll}\exp[-(v_{t}-v_{0}+\mu
t)^{2}/(4Dt)-v_{t}/D]&&\mbox{for }v_{t}<v^{-}(t)\\\ \exp[-(v_{t}-v_{0}+\mu
t)^{2}/(4Dt)]&&\mbox{for }v_{t}>v^{+}(t)\\\ \exp(-\mu|v_{t}|/D)&&\mbox{for
}v_{t}\in[v^{-}(t),v^{+}(t)],\end{array}\right.$ (21)

where $N_{1}$ and $N_{2}$ are normalization factors. This result is compared
in Fig. 3 with the exact propagator. There we see that the SPA(0) is a good
approximation of $p(v,t|v_{0},0)$ at short and long times, but does not
capture the bimodal structure of the propagator arising when the optimal path
hits the origin for times close to $t=v_{0}/\mu$. While a kink of the exact
propagator shows up at time $t<v_{0}/\mu$ in the region $v_{t}>0$ [see Fig.
3(b)], the corresponding kink of the SPA(0) only appears at time $t>v_{0}/\mu$
[see Fig. 3(c) as well as Eqs. (20) and (21)]. For comparison, we also show in
Fig. 3 the result of the propagator obtained from Langevin simulations of the
SDE using the standard Euler-Maruyama integration scheme. The application of
this scheme is stable for the piecewise-smooth SDE and only requires that we
choose the integration time step small enough so that we can reproduce the
cusp seen in its propagator.

Figure 3: (Color online) Propagator of the dry friction model (1) for initial
velocity $v_{0}=2$, $\mu=1$, $D=0.01$ and different values of time: (a) $t=1$,
(b) $t=1.8$, (c) $t=2.2$, and (d) $t=3$. Here DF denotes the exact analytical
result (3), MCS the Monte Carlo simulation of Eq. (1) with step size $0.0001$
and an ensemble of $10^{6}$ realizations, SPA(0) the leading-order saddle-node
expansion [see Eqs. (20) and (21) in Sec. II.1)], and SPA(1) the higher-order
saddle-node expansion [see Eqs. (25), (27), and (29) in Sec. II.2)].

We will see in the next section that the inclusion of the Jacobian in the
saddle-point approximation enables us to reproduce this bimodality more
accurately. An important remark before we get to that part is that the SPA(0)
yields the same approximation as that obtained from the Freidlin-Wentzell (FW)
large-deviation theory of dynamical systems perturbed by noise Freidlin and
Wentzell (1984). Indeed, it is not difficult to verify that the limit

$I(v,t|v_{0},0)=\lim_{D\rightarrow 0}-4D\,\ln
p(v,t|v_{0},0)=\lim_{D\rightarrow
0}-4D\,\ln\hat{p}\left(\frac{\mu}{D}v_{t},\frac{\mu^{2}}{D}t\,\right|\\!\left.\frac{\mu}{D}v_{0},0\right),$
(22)

which defines in the FW theory the so-called pseudo- or quasipotential
$I(v,t|v_{0},0)$, is equal to the SPA(0) action. This shows that the
application of the FW theory to nonsmooth SDEs can be done following the
heuristic principles formulated above. In this case, the different optimal
paths that we have found can be associated, in the low-noise limit, with two
different physical modes of motion of the noiseless system: Direct paths
represent slip motion, whereas indirect paths represent stick motion Baule et
al. (2010, 2011). For any fixed value of $v_{t}$, the optimal path will always
be an indirect path if the time is sufficiently large, which means that all
optimal paths are indirect paths in the limit $t\to\infty$.

###  II.2 Corrected action

One way to correct the SPA(0) is to use the optimal paths obtained before but
to evaluate their action by including the Jacobian in the action:

$S[v]=\int_{0}^{t}\left\\{[\dot{v}+\mu\sigma(v)]^{2}-4D\mu\delta(v)\right\\}\,\mathrm{d}s.$
(23)

This defines a first-order saddle-point approximation [SPA(1)] of the
propagator that retains the optimal paths of SPA(0) but includes the
subdominant correction of the Jacobian, which is multiplied by $D$.

To obtain this approximation, we need to evaluate how the contribution arising
from the $\delta$ function in Eq. (23) changes the action for the different
paths considered. For direct paths, there is obviously only a contribution
when the $v=0$ axis is crossed, so that

$\int_{0}^{t}\delta(v_{d}^{(0)})\mathrm{d}s=\left\\{\begin{array}[]{lll}t/(v_{0}-v_{t})&&\mbox{for
}v_{t}<0\\\ 0&&\mbox{for }v_{t}>0.\\\ \end{array}\right.$ (24)

Thus, the corresponding corrected action is

$S[v_{d}^{(0)}]=\left\\{\begin{array}[]{lcl}(v_{t}-v_{0}+\mu t)^{2}/t-4\mu
v_{t}-4D\mu t/(v_{0}-v_{t})&\mbox{ for }&v_{t}<0\\\ (v_{t}-v_{0}+\mu
t)^{2}/t&\mbox{ for }&v_{t}>0.\end{array}\right.$ (25)

For intermediate paths, the evaluation of the Jacobian term is straightforward
as well. Any sensible representation of the $\delta$ function will result in a
symmetric average of the inverse slopes of the path, leading to

$\int_{0}^{t}\delta(v_{int}^{(0)})ds=\frac{1}{2\big{|}\dot{v}^{(0)}_{int}(t_{0}^{-},t_{0})\big{|}}+\frac{1}{2\big{|}\dot{v}^{(0)}_{int}(t_{0}^{+},t_{0})\big{|}}=\frac{t}{v_{0}+|v_{t}|},$
(26)

where we have used the condition $t_{0}=v_{0}t/(v_{0}+|v_{t}|)$.  Hence it
follows from Eqs. (23) and (26) that

$S[v_{int}^{(0)}]=(|v_{t}|+v_{0}-\mu t)^{2}/t+4\mu|v_{t}|-4D\mu
t/(v_{0}+|v_{t}|).$ (27)

Indirect paths require a closer inspection: These paths vanish over an entire
interval of time, so the contribution originating from the Jacobian is ill
defined. To treat this, we follow the previous principle that paths on the
attractor do not contribute to the action and posit the following two
additional heuristic principles: (iii) Parts of indirect paths on the
attractor are not considered as contributing to the Jacobian and (iv)
nonvanishing parts of indirect paths contribute, as for intermediate paths, to
the Jacobian in a weighted average way. With these principles, the corrected
action of indirect paths is finite and equal to

$\int_{0}^{t}\delta(v_{ind}^{(0)})\mbox{d}s=\frac{1}{2\big{|}\dot{v}_{ind}^{(0)}(t_{1}^{-},t_{1},t_{2})\big{|}}+\frac{1}{2\big{|}\dot{v}_{ind}^{(0)}(t_{2}^{+},t_{1},t_{2})\big{|}}=\mu,$
(28)

where we have used the conditions $t_{1}=v_{0}/t$ and $t_{2}=t-|v_{t}|/\mu$.
Thus the corresponding action (23) can be evaluated as

$S[v_{ind}^{(0)}]=4\mu|v_{t}|-4D\mu\qquad\mbox{ for }t>(v_{0}+|v_{t}|)/\mu.\\\
$ (29)

By properly comparing the corrected actions (25), (27), and (29), we can
determine which path is minimal depending on $v_{t}$ and $t$ to obtain the
propagator $p^{(0;1)}(v,t|v_{0},0)$ [see Eq. (75)]. The result of this
minimization is shown in Fig. 3 as SPA(1). We see that the inclusion of the
Jacobian correction qualitatively improves the propagator as compared to the
lowest-order approximation, SPA(0) of Sec. II.1, even though there are still
some deviations in the transient regime where the exact propagator shows a
bimodality. In particular, at time $t<v_{0}/\mu$ the SPA(1) already shows the
kink corresponding to that of the exact propagator in the region $v_{t}>0$
[see Fig. 3(b)]. This phenomenon is different with that of the SPA(0), which
does not have such a kink at $t<v_{0}/\mu$. We also observe in Fig. 3(c) that
the SPA(1) develops two kinks in the region $v_{t}>0$ in an intermediate time
interval. These kinks are caused by the intermediate path, which becomes the
optimal path for certain choices of the end points. The SPA(0), by contrast,
has only one kink in the positive region, since for this approximation
intermediate paths can never be the optimal path [see Eq. (18)]. These kinks
are a well-known artifact of low-noise approximations and are smoothed in the
exact propagator by the finite diffusion.

###  II.3 Corrected action with corrected path

The SPA(1) corrects the SPA(0) by including the Jacobian term in the action,
while using the optimal paths of SPA(0), i.e., the paths that minimize the
zeroth-order action $S^{(0)}[v]$. As a further correction to this
approximation, it is tempting to obtain the optimal paths by minimizing the
corrected action $S[v]$ with the Jacobian, thereby constructing a “full”
first-order approximation.

Unfortunately, this approach does not work as the action turns out to diverge
for $v_{t}=0$. To see this, evaluate the action of Eq. (23) for an
intermediate path $v_{int}^{(0)}$ with a kink at $t_{0}=t/2$ (see Fig. 1):

$S[v_{int}^{(0)}]=2\mu(|v_{t}|-v_{0})+\mu^{2}t+2(v_{0}^{2}+v_{t}^{2})/t-D\mu\left(t/v_{0}+t/|v_{t}|\right).$
(30)

This value is an upper bound for the minimum of the action, which determines
the density according to Eq. (77). The problem with this result is that
$S[v_{int}^{(0)}]\rightarrow-\infty$ when $v_{t}\rightarrow 0$, leading to a
(nonintegrable) singularity for the propagator $p(v_{t},t|v_{0},0)$ at
$v_{t}=0$. Thus the approximation scheme based on obtaining the optimal paths
from $S[v]$ results in a non-normalizable expression, which implies that a
full first-order SPA is not possible for the SDE of Eq. (1). It is clear that
this problem will also arise in any SDE having, as in Eq. (1), points where
the force of the SDE is discontinuous. One way to approach this problem is to
explore regularizations of such discontinuities.

##  III Regularized SDE

As seen in the previous section, the weak-noise expansion of the path integral
for nonsmooth SDEs faces some difficulties related to the minimization of the
action and the singularity of the Jacobian term. To treat this problem, we now
consider the regularized SDE of Eq. (7) in which the discontinuous drift
$\sigma(v)$ is replaced by the smooth drift $\tanh\left(v/\varepsilon\right)$
involving the additional (small) parameter $\varepsilon$. For this smooth SDE,
the aforementioned difficulties do not occur: We can minimize the action over
smooth differentiable paths and the Jacobian is well defined. In this section,
we are interested in understanding the weak-noise properties of this
regularized SDE.

To investigate the regularized model, we introduce non-dimensional units:

$u=v/\varepsilon,\qquad\tau=\mu t/\varepsilon.$ (31)

Equation (7) is then simply written as

$\dot{u}=-\tanh(u)+\sqrt{\widetilde{D}}\,\xi(\tau),$ (32)

where $\widetilde{D}=D/(\varepsilon\mu)$ is now the only parameter of the
model, called the effective diffusion constant. An important point to note is
that the two limits $\varepsilon\rightarrow 0$ and $D\rightarrow 0$ do not
commute. In the following, we will be interested in the smooth model for small
and moderate effective noise amplitudes.

As mentioned in the Introduction, the propagator of the regularized model does
not have a known closed analytic form. Thus, for benchmarking the weak noise
expansion results, we need to resort to numerical simulations obtained with
the Euler-Maruyama scheme Kloeden and Platen (1992), which accurately
reproduce, as for the piecewise smooth SDE, all the features of the
propagator, provided we choose the integration time step small enough.

As before, the leading order of the weak-noise expansion of Eq. (32) is
determined by the action [see Eqs. (72 and (74)]

$S^{(0)}[u]=\int_{0}^{\tau}\left[\dot{u}(s)+\tanh u(s)\right]^{2}\mbox{d}s$
(33)

evaluated for the solution of the Euler-Lagrange boundary value problem [see
Eq. (73)]

$\ddot{u}^{(0)}(s)=\tanh[u^{(0)}(s)]/\cosh^{2}[u^{(0)}(s)],\quad
u^{(0)}(0)=u_{0},\quad u^{(0)}(\tau)=u_{\tau}.$ (34)

The leading order may be improved by consistently taking first-order
contributions into account. The propagator (77) is then determined by the
action [see Eq. (71)]

$S[u]=\int_{0}^{\tau}\left\\{\left[\dot{u}(s)+\tanh
u(s)\right]^{2}-2\widetilde{D}/\cosh^{2}[u(s)]\right\\}\mbox{d}s$ (35)

evaluated at the first-order path that obeys the boundary value problem (76)

$\ddot{u}^{(1)}(s)=\big{(}1+2\widetilde{D}\big{)}\tanh[u^{(1)}(s)]/\cosh^{2}[u^{(1)}(s)],\quad
u^{(1)}(0)=u_{0},\quad u^{(1)}(\tau)=u_{\tau}.$ (36)

The two boundary value problems (34) and (36) differ just by a rescaling of
time and can be solved using a numerical shooting method (see Appendix A).

Figure 4 compares the numerical Monte Carlo simulations of Eq. (32) with the
zeroth- and first-order weak-noise approximations for small and moderate noise
amplitudes. We see that, while the short- and long-time behaviors are captured
quite well by the lowest order approximation, substantial transient deviations
are visible when the maximum of the propagator approaches the origin. The
first-order approximation is in fact able to deal with such a feature, even
for substantial noise amplitudes. Thus the scheme outlined above can be
considered as a candidate to deal with the weak-noise limit even in systems
that mimic the discontinuous drift.

Figure 4: (Color online) Propagator of Eq. (32) for small and moderate values
of the effective diffusion, (a)–(c) $\widetilde{D}=0.01$ and (d)–(f)
$\widetilde{D}=0.1$, initial condition $u_{0}=10$ and three different values
of $\tau$. Here MCS denotes the Monte Carlo simulations of the Langevin
equation (32) with time step $0.005$ and an ensemble of $10^{6}$ realizations,
SPA(0) the lowest order of the saddle node expansion using the action of Eq.
(33) [see Eq. (74)], and SPA(1) the first order of the saddle-node expansion
using the action of Eq. (35) [see also Eq. (77)].

Naturally, from our study of the discontinuous model, we cannot expect the
weak-noise approximation to yield the full propagator of the smooth model in
the asymptotic limit $\varepsilon\rightarrow 0$. However, it should be
possible to obtain, at a quantitative level, the main features of the
propagator for suitable small values of $D$ and $\varepsilon$. Of course, an
improved scheme such as the first-order expansion needs to be applied, as one
cannot rely on extremely small values of the effective diffusion to cover
cases that are sufficiently close to a discontinuous drift. Indeed, Fig. 5
shows that this is the case if we translate the results obtained via Eqs. (35)
and (36) to the dimensional units via Eq. (31). The SPA(1) performs well in
short time and is also able to capture the main profiles of the piecewise-
smooth SDE in moderate and long times. Larger deviations between the
discontinuous and the regularized result appear only in a neighborhood of size
$\varepsilon$ of the discontinuity. Hence we may conclude that, in these
particular cases, a suitable regularization and first-order saddle-point
approximation are able to capture the essential features of the piecewise
smooth SDE.

Figure 5: (Color online) Comparison of the propagators of the dry friction
model (1) and of the regularized scheme (7) for $\mu=1$, $D=0.01$,
$\varepsilon=0.05$, initial condition $v_{0}=2$, and three different values of
time: (a) $t=1$, (b) $t=1.9$, and (c) $t=3$. Here DF denotes Analytical result
for the propagator of the dry friction model [see Eq. (3)], MCS the Monte
Carlo simulation of the Langevin equation (7) with time step $0.005$ and an
ensemble of $10^{6}$ realizations, SPA(1) the first order saddle-point
approximation scheme using the action of Eq. (35). The inset in (b) shows the
SPA(1) of the smooth case (solid line) and the SPA(1) of the dry friction case
(dashed line) [see Eqs. (25), (27), and (29)]. Both of these graphs have a
kink that matches the structure appearing in the exact solution of the dry
friction model (DF). The inset in (c) shows the data for the DF and SPA(1) on
a larger scale on a semilogarithmic plot. Large quantitative deviations appear
only in a neighborhood of size $\varepsilon$ of the discontinuity.

To close this section, let us comment on the kinks appearing in the
quasipotential, corresponding to the minimized action. These nonanalytic
points are well known in large deviation theory to arise from the
nonuniqueness of the solutions of the Euler-Lagrange boundary value problem
Graham and Tél (1985). Since Eqs. (34) and (36) differ only by a rescaling of
time, and since we are here mainly concerned with the essential structure of
the quasipotential we just focus on the simple zeroth-order expansion (33) and
(34). Figure 6 shows the solution of the boundary value problem using the
shooting method for a given value of $u_{0}$. At some finite value of $\tau$,
the boundary value problem develops a cusp singularity, beyond which three
solutions occur in a finite interval of $u_{\tau}$ values, corresponding to
smooth versions of the previously identified direct, indirect, and
intermediate paths. Hence we recover in the smooth model the same path
structure of the discontinuous model.

Figure 6: (Color online) Solution of the boundary value problem (34) with the
shooting method for the given value $u_{0}=4$ and three different values of
$\tau$. (a)–(c) Solutions of the differential equation (34) for different
values of the initial slope $\dot{u}^{(0)}(0)=p_{u}(0)$. (d)–(f) Dependence of
the final value $u_{\tau}=u^{(0)}(\tau)$ on the initial condition
$\dot{u}^{(0)}(0)=p_{u}(0)$ of the differential equation (34), indicated by
the solid line. The symbols indicate the solution that minimizes the action
(33).

A nonunique solution to the boundary value problem implies that the actual
optimizing path has to be determined from minimizing Eq. (33) among the three
possible paths. In our case, Fig. 6 shows that the optimal path is either a
direct or an indirect path, as for the piecewise-constant SDE. The
intermediate path is always a saddle point of the action, which implies that
the quasipotential consists of two analytic branches, with a kink appearing
when the type of optimal paths changes [see Figs. 5(b), 7(b), and 7(c)]. The
kink, which also shows up in the numerical solution of the Fokker-Planck
equation, is an important feature of the quasipotential: It appears at finite
time and then moves to larger $u_{\tau}$ values [see Fig. 7(c)]. This feature,
which is related to the convergence of the propagator towards the stationary
distribution, is studied in detail next.

Figure 7: (Color online) quasipotential of the regularized model, Eq. (32),
for initial condition $u_{0}=4$, $D=0.0001$ and three different values of
time: (a) $\tau=4$, (b) $\tau=7$, and (c) $\tau=10$. Here FP (dashed line)
denotes the potential $-4D\ln p(u_{\tau},\tau|u_{0},0)$ computed by numerical
integration of the Fokker-Planck equation and SPA(0) (solid line) the result
of the leading-order SPA, i.e., minimized action (33). The inset in (c) shows
the position of the kink, $u_{k}$, as a function of the time $\tau$, as
obtained from the SPA(0). The closed circle shows the time and the position
where the kink emerges [see also Fig. 12].

##  IV Analytic properties of the action

The previous analysis was mainly based on a numerical solution of the boundary
value problem (34) and the evaluation of the corresponding action integral
(33) (see Cassing and Knoll (1981) for a related numerical study). In this
section, we obtain further insights into the weak-noise limit by studying some
analytical properties of the quasipotential, focusing on the zeroth-order
approximation of the action. The first-order approximation of the action can
be analyzed along similar lines.

###  IV.1 Action integral

Figure 8: Phase portrait of the Hamiltonian system (37) that governs the
boundary value problem (34). The region of positive energy is shaded.

The Euler-Lagrange equations (34) are equivalent to a conservative system with
the Hamiltonian

$H(p_{u},u)=\frac{1}{2}p_{u}^{2}-\frac{1}{2}\tanh^{2}(u),$ (37)

where $p_{u}=\dot{u}$. The phase portrait of this system, shown in Fig. 8, is
important to understand the structure of the boundary value problem. The
origin in phase space is a hyperbolic equilibrium point. Regions of positive
and negative energy, respectively, are bounded by the separatrix of this fixed
point. Solutions of the equation of motion (34) are given by constant energy
levels

$\dot{u}=p_{u}=\pm\sqrt{\tanh^{2}(u)+2H},$ (38)

where the sign determines whether the path is in the positive or negative
momentum region of phase space. The zeroth-order action (33) for a path with
energy $H$, starting at $u_{0}$ and terminating at $u_{\tau}$ can then be
written as

$S^{(0)}[u]=S(u_{0},u_{\tau},\tau,H)=2\int_{u_{0},H}^{u_{\tau}}p_{u}(s)\mbox{d}u(s)+2\ln\left(\frac{\cosh(u_{\tau})}{\cosh(u_{0})}\right)-2H\tau$
(39)

if we take the expression (37) for the Hamiltonian into account. The duration
$\tau$ of this path from $u_{0}$ to $u_{\tau}$ is evaluated as

$\tau=\int_{u_{0},H}^{u_{\tau}}\frac{\mbox{d}u(s)}{p_{u}(s)},$ (40)

and allows us to express the energy in Eq. (39) in terms of the time. Thus the
problem has been reduced to evaluating two integrals and solving algebraic
equations. To take care of the correct sign in Eq. (38), we have to
distinguish three cases, according to the sign of the terminal point
$u_{\tau}$.

####  IV.1.1 Case 1: $u_{\tau}<0<u_{0}$

As can be seen from the phase portrait of Fig. 9(a), there is for given energy
$H>0$ a unique path connecting the two boundary points.

Figure 9: (Color online) (a) Phase portrait of the Hamiltonian (37) with an
orbit connecting two boundary points $u_{\tau}<0<u_{0}$ at a given energy
$H>0$. (b) Energy time relation (41) for $u_{0}=10$ and $u_{\tau}=-2$.

Along this path Eq. (38) holds with a minus sign, and the expression relating
energy and time, Eq. (40), results in

$\tau=\int_{u_{\tau}}^{u_{0}}\frac{\mbox{d}u}{\sqrt{\tanh^{2}(u)+2H}}=\theta_{-}(u_{0},u_{\tau},H),$
(41)

where

$\theta_{-}(u_{0},u_{\tau},H)=\frac{1}{\sqrt{1+2H}}\ln\left(\frac{\sinh(u_{0})+\sqrt{2H/(1+2H)+\sinh^{2}(u_{0})}}{\sinh(u_{\tau})+\sqrt{2H/(1+2H)+\sinh^{2}(u_{\tau})}}\right).$
(42)

It is easy to show that in the range $H>0$ and $u_{\tau}<0<u_{0}$, Eq. (42) is
a monotonically decreasing function of $H$ [see Fig. 9(b)], i.e., Eq. (41)
defines the energy $H$ as an analytic expression of $\tau$ and of the boundary
points (see Appendix B.1). As for the integral that enters the action (39) we
obtain

$\int_{u_{0},H}^{u_{\tau}}p_{u}(s)du(s)=\int_{u_{\tau}}^{u_{0}}\sqrt{\tanh^{2}(u)+2H}\mbox{d}u=\int_{u_{\tau}}^{u_{0}}\frac{1+2H-1/\cosh^{2}(u)}{\sqrt{\tanh^{2}(u)+2H}}\mbox{d}u=(1+2H)\tau-\sigma_{-}(u_{0},u_{\tau},H),$
(43)

where we have introduced the abbreviation

$\sigma_{-}(u_{0},u_{\tau},H)=\ln\left(\tanh(u_{0})+\sqrt{2H+\tanh^{2}(u_{0})}\right)-\ln\left(\tanh(u_{\tau})+\sqrt{2H+\tanh^{2}(u_{\tau})}\right).$
(44)

Therefore, the entire action (39) finally reads

$S(u_{0},u_{\tau},\tau,H)=2(1+H)\tau+2\ln\left(\frac{\cosh(u_{\tau})}{\cosh(u_{0})}\right)-2\sigma_{-}(u_{0},u_{\tau},H).$
(45)

Since the boundary value problem has a unique solution, the inverse of the
energy-time relation (41) is single valued, so that the expression for the
action defines an analytic expression of time $\tau$ and of the boundary
points in the region $u_{\tau}<0<u_{0}$ (see as well Fig. 7).

Figure 10: (Color online) (a) Phase portrait of the Hamiltonian (37) with
direct (1) and indirect (2) paths connecting boundary points
$0<u_{\tau}<u_{0}$. (b) Energy time relation (41) (solid line) and (46)
(dashed line) for small values of $u_{0}$ ($u_{0}=2$ and $u_{\tau}=1$) with a
unique inverse. (c) Energy time relations (41) (solid line) and (46) (dashed
line) for larger values of $u_{0}$ ($u_{0}=3$ and $u_{\tau}=2$) with a
multivalued inverse in a finite time interval.

####  IV.1.2 Case 2: $0<u_{\tau}<u_{0}$

For boundary points $0<u_{\tau}<u_{0}$ the phase portrait in Fig. 10(a) and
Eq. (37) show that the connecting path has energy
$H>H_{min}=-\tanh^{2}(u_{\tau})/2$. On the one hand, for energies $H>H_{min}$,
there exists a direct path whose duration is given by Eq. (41) and whose
action is determined by Eq. (45). On the other hand, for $H_{min}<H<0$, there
exists an indirect path with a turning point at
$u_{min}=\mbox{artanh}(\sqrt{-2H})$. The duration of this path is obtained
from Eq. (40) by splitting the path into two parts and using appropriate signs
in Eq. (38):

$\displaystyle\tau=\int_{u_{0}}^{u_{min}}\frac{\mbox{d}u}{-\sqrt{\tanh^{2}u+2H}}+\int_{u_{min}}^{u_{\tau}}\frac{\mathrm{d}u}{\sqrt{\tanh^{2}u+2H}}=\theta_{+}(u_{0},u_{\tau},H),$
(46)

where

$\displaystyle\theta_{+}(u_{0},u_{\tau},H)$ $\displaystyle=$
$\displaystyle\frac{\ln\left(\sinh(u_{0})+\sqrt{2H/(1+2H)+\sinh^{2}(u_{0})}\right)}{\sqrt{1+2H}}+\frac{\ln\left(\sinh(u_{\tau})+\sqrt{2H/(1+2H)+\sinh^{2}(u_{\tau})}\right)}{\sqrt{1+2H}}$
(47) $\displaystyle+\frac{\ln\left((1+2H)/(-2H)\right)}{\sqrt{1+2H}}.$

Similarly, the integral appearing in the action (39) can be evaluated [see Eq.
(43)]

$\displaystyle\int_{u_{0},H}^{u_{\tau}}p_{u}(s)\mbox{d}u(s)=-\int_{u_{0}}^{u_{min}}\sqrt{2H+\tanh^{2}(u)}\mbox{d}u+\int_{u_{min}}^{u_{\tau}}\sqrt{2H+\tanh^{2}(u)}\mbox{d}u=(1+2H)\tau-\sigma_{+}(u_{0},u_{\tau},H),$
(48)

where we have introduced the abbreviation

$\displaystyle\sigma_{+}(u_{0},u_{\tau},H)=\ln\left(\tanh(u_{0})+\sqrt{2H+\tanh^{2}(u_{0})}\right)+\ln\left(\tanh(u_{\tau})+\sqrt{2H+\tanh^{2}(u_{\tau})}\right)-\ln(-2H).$
(49)

Thus, for the action (39) of this path we obtain

$S(u_{0},u_{\tau},\tau,H)=2(1+H)\tau+2\ln\left(\frac{\cosh(u_{\tau})}{\cosh(u_{0})}\right)-2\sigma_{+}(u_{0},u_{\tau},H).$
(50)

To find the minimizing path, we have to take a closer look at the energy-time
relation. This relation consists of two branches [see Fig. 10(b) and 10(c)].
The branch determined by Eq. (41) is defined for all energies $H_{min}<H$ and
is a monotonically decreasing function for $0<u_{\tau}<u_{0}$ (see Appendix
B.1). The second branch, defined for negative energies $H_{min}<H<0$ only, is
given by Eq. (46). This branch is monotonically increasing for small values of
$u_{0}$ and develops an inflection point if the boundary points exceed a
critical value, which implies that, for small values of $u_{0}$, the relation
determines the energy in terms of time and boundary points uniquely. This
uniqueness is expected, as the analysis of the regularized model reduces, in
the limit of small value of $u_{0}$, to that of the Ornstein-Uhlenbeck
process. However, if the initial condition is beyond the linear regime, the
energy time relation has a multivalued inverse within a finite-time interval
[see Fig. 10(c)]. The value of the action entering in the propagator is then
the minimum among the three possible inverse values. The data show that for
fixed values of $u_{0}$ and $u_{\tau}$ the minimum is given by the low energy
solution, i.e., by the direct path, up to a critical value of $\tau$, while
the minimum for larger values of $\tau$ is given by the high-energy solution,
i.e., the indirect path. The intermediate path, i.e., the solution with energy
in between, is always of saddle type. As a consequence, the action, considered
as a function of $u_{\tau}$ switches the analytical branch and develops a kink
[see Figs. 7(b), 7(c), and 11(b)].

Figure 11: (Color online) (a) Phase portrait of the Hamiltonian (37) with
direct (1) and indirect (2) orbits connecting boundary points
$0<u_{0}<u_{\tau}$. (b) Minimal action (39) as a function of $u_{\tau}$ in the
interval $u_{\tau}>u_{0}=10$ for different values of the time $\tau$, showing
a kink moving to the right.

####  IV.1.3 Case 3: $0<u_{0}<u_{\tau}$

As in the previous case there are two types of paths to consider: In the
energy range $H>H_{min}=-\tanh^{2}(u_{0})/2$, there exists a direct path
connecting the two boundary points, while for negative energies $H_{min}<H<0$,
there exist in addition indirect paths with a turning point at
$u_{min}=\mbox{artanh}(\sqrt{-2H})$. This is similar to case 2 and the
relevant integrals can be dealt with by applying the symmetry of the phase
portrait [see Fig. 11(a)]. To find the duration of the direct path (40) we
observe that just by swapping the boundary points and reversing the direction
of the path we obtain a corresponding mirror orbit with the same duration,
which we have dealt with in the previous paragraph. Thus Eq. (41) tells us
that

$\tau=\int_{u_{0}}^{u_{\tau}}\frac{\mbox{d}u}{\sqrt{\tanh^{2}u+2H}}=\theta_{-}(u_{\tau},u_{0},H).$
(51)

Since the integral of the action has the same symmetry, we obtain the result
stated in Eq. (43) by swapping $u_{0}$ and $u_{\tau}$ on the right-hand side.
Thus the action (39) finally reads

$S(u_{0},u_{\tau},\tau,H)=2(1+H)\tau+2\ln\left(\frac{\cosh(u_{\tau})}{\cosh(u_{0})}\right)-2\sigma_{-}(u_{\tau},u_{0},H).$
(52)

The action itself does not share the aforementioned symmetry because of the
additional terms appearing in Eq. (39).

For the indirect path the same reasoning applies. The duration of the path,
given by Eq. (46), applies since the right-hand side is symmetric in $u_{0}$
and $u_{\tau}$ and swapping both arguments does not have any effect. In fact,
the same argument is valid for the integral that enters the action. Therefore,
Eq. (48) holds since the right-hand side is a symmetric expression in the
boundary points and the corresponding action is indeed given by Eq. (50).

Overall, we see that the energy-time relationship is at the heart of
understanding the structure of the optimal paths, the minimized action, and
finally the propagator of the SDE. As in the previous case, the relation
consists of two branches [see Figs. 10(b) and 10(c)]. The lower branch is
given by the direct path (51), which is a monotonically decreasing function on
its domain $H>H_{min}$. The second branch for negative energies $H_{min}<H<0$
is determined by Eq. (46) and thus largely the discussion of the previous
section applies. If the initial condition is small, so that no cubic
singularity appears for $0<u_{\tau}<u_{0}$ (case 2), then such a singularity
will finally develop for sufficiently large value of $u_{\tau}$ resulting in a
non analytic minimized action. If, on the contrary, $u_{0}$ is so large that
the occurrence of the kink has already happened in the domain
$0<u_{\tau}<u_{0}$, then the kink just propagates to larger values of
$u_{\tau}$ [see Fig. 11(b)].

###  IV.2 Asymptotics of the action

The analytic expressions derived in the previous section allow us to study in
some detail the properties of the propagator shown, e.g., in Fig. 7. Of
particular relevance is the approach to the stationary state, the tail
behavior of the distribution, and the emergence and dynamics of
nonanalyticities of the quasipotential.

####  IV.2.1 Approach to equilibrium

To study the convergence of the propagator towards the stationary
distribution, we consider a fixed value of $u_{0}$ and $u_{\tau}$ and take the
asymptotic limit $\tau\rightarrow\infty$. It is obvious from the phase
portrait or from the energy-time relation [see Figs. 9(b) and 10(b)] that the
asymptotic limit implies $H\rightarrow 0$. Hence the energy-time relation is
determined by Eq. (41) or (46), depending on the sign of $u_{\tau}$. In both
cases, a straightforward expansion yields

$\tau=\ln(2\sinh(u_{0}))+\ln(-\sinh(u_{\tau})/H)+O(H\ln|H|)\,.$ (53)

For the action, either Eq. (45) or (50) applies, depending on the sign of
$u_{\tau}$, yielding the expansion

$\displaystyle
S(u_{0},u_{\tau},\tau,H)=2(1+H)\tau+2\ln\left(\frac{\cosh(u_{\tau})}{\cosh(u_{0})}\right)-2\left(\ln\left(2\tanh(u_{0})\right)+\ln\left(-\tanh(u_{\tau})/H\right)\right)+O(H).$
(54)

By solving Eq. (53) to leading order for $H$, Eq. (54) yields, as expected,
the potential of the stationary distribution [see also Eq. (74)]

$S(u_{0},u_{\tau},\tau,H)=4\ln\cosh(u_{\tau})+O(\tau\exp(-\tau)).$ (55)

The full action as well as the propagator at finite time depends on the
initial condition $u_{0}$ and hence does not have a symmetry with respect to
$u_{\tau}$. However, the stationary action and the corresponding stationary
distribution is independent of $u_{0}$ [see Eq. (55)] in the limit
$\tau\rightarrow\infty$. In that limit the symmetry under the change
$u_{\tau}\rightarrow-u_{\tau}$ is restored [see as well Eq. (6)]. Above all,
the saddle-point expansion preserves all the symmetries of the underlying
equations of motion. As for the transient properties it is easy, though more
tedious, to compute the subleading term in Eq. (55), which then quantitatively
describes the relaxation process.

####  IV.2.2 Tail behavior of the propagator

To investigate the action for large values of arguments, we consider the
asymptotic limit $|u_{\tau}|\rightarrow\infty$ for fixed values of the time
$\tau$ and initial condition $u_{0}$, which implies, from the phase portraits
of Figs. 9(a) and 11(a), $H\rightarrow+\infty$. The energy-time relation is
determined by Eqs. (41) and (51), depending on the sign of $u_{\tau}$, and a
direct expansion results in the expression

$\tau=\frac{1}{\sqrt{1+2H}}\left(|u_{\tau}-u_{0}|+O(H^{-1})\right),$ (56)

which is valid irrespective of the sign of $u_{\tau}$. This expression is
easily inverted to obtain

$1+2H=(u_{\tau}-u_{0})^{2}/\tau^{2}+O(u_{\tau}^{-1}).$ (57)

For the action Eqs. (45) or (52) apply, depending on the sign of $u_{\tau}$.
Expansion in terms of $H$ yields $\sigma_{-}$ to be of order $O(H^{-1})$ so
that

$S(u_{0},u_{\tau},\tau,H)=2(1+H)\tau+2\ln\left(\frac{\cosh(u_{\tau})}{\cosh(u_{0})}\right)+O(H^{-1}).$
(58)

Hence, using Eq. (57) we end up with

$S(u_{0},u_{\tau},\tau,H)=\frac{(u_{\tau}-u_{0})^{2}}{\tau}+\tau+2\ln\left(\frac{\cosh(u_{\tau})}{\cosh(u_{0})}\right)+O(u_{\tau}^{-1}).$
(59)

To leading order, the tails of the propagator have a Gaussian shape, as
expected, while the subleading terms indicate a transition from the Gaussian
shape to the exponential shape of the stationary distribution.

####  IV.2.3 Emergence and propagation of the kink

As mentioned before, a kink in the minimized action appears when the energy-
time relation has multivalued inverse. Since
$\theta_{-}(u_{0},u_{\tau},\tau,H)$ is a monotonically decreasing function of
$H$, the kink appears when $\theta_{+}(u_{0},u_{\tau},\tau,H)$ develops a
crossover from a monotonically increasing shape into a cubic shape [see Figs.
10(b) and 10(c)]. Hence the time $\tau_{c}$ when the kink appears and its
position $u_{\tau}$ as a function of initial position $u_{0}$ are determined
by

$\tau=\theta_{+}(u_{0},u_{\tau},H),\qquad\frac{\partial\theta_{+}}{\partial
H}=0,\qquad\frac{\partial^{2}\theta_{+}}{\partial H^{2}}=0.$ (60)

Figure 12 shows the solution of these equations. For small values of $u_{0}$,
i.e., in the limiting case of the Ornstein-Uhlenbeck process, both the time
and the position for the emergence of the kink become large. For large values
of $u_{0}$, i.e., when the drift is almost constant, the time becomes large
again but the kink appears close to the origin. At an intermediate scale of
the size of the regularization length the time for the occurrence of the kink
becomes minimal and the kink appears close to the initial condition.

Figure 12: (Color online) Time $\tau_{c}$ and position $u_{\tau}$ for the
emergence of the kink in the minimized action as a function of the initial
condition $u_{0}$ [see Eqs. (60)]. The closed circle indicates the starting
point of the kink for $u_{0}=4$, as in Fig. 7(c).

After its appearance the kink starts to move with positive velocity [see Figs.
11(b) and 7(c)]. We can quantify such a dynamical feature in terms of an
asymptotic expansion for large $\tau$ and large $u_{\tau}$, keeping the ratio
$c=u_{\tau}/\tau$ of $O(1)$. The key to the analysis is the shape of the
energy-time relation in such a limit (see Fig. 13).

Figure 13: (Color online) Energy time relation (51) and (46) for $u_{0}=2$ and
$u_{\tau}=5$. The action takes its minimum on the lower branch $AB$ if
$\tau<10.5$ and on the upper branch $CD$ if $\tau>10.5$.

The relation consists of the monotonic lower branch (51) determined by
$\theta_{-}$ and a cubic shaped upper branch (46) determined by $\theta_{+}$.
Depending on the value of $\tau$ the action takes its minimum either on the
lower branch at an energy $H_{-}$ that stays of order $O(1)$ or on the upper
branch at an energy $H_{+}$ that tends to zero in the asymptotic limit
considered here. Hence Eqs. (51) and (46) yield

$\displaystyle\tau$ $\displaystyle=$
$\displaystyle\theta_{-}(u_{\tau},u_{0},H_{-})=\frac{u_{\tau}}{\sqrt{1+2H_{-}}}+O(1),$
(61) $\displaystyle\tau$ $\displaystyle=$
$\displaystyle\theta_{+}(u_{\tau},u_{0},H_{+})=u_{\tau}+\ln(2\sinh(u_{0}))-\ln(-2H_{+})+O(\exp(-2u_{\tau}),H_{+}).$
(62)

Both expressions are easily inverted to obtain

$1+2H_{-}=(u_{\tau}/\tau)^{2}+O(\tau^{-1})$ (63)

and

$-2H_{+}=2\sinh(u_{0})\exp(u_{\tau}-\tau)\left(1+O(\exp(-2u_{\tau}),\exp(-\tau(1-c)))\right),\quad(c=u_{\tau}/\tau<1),$
(64)

where for Eq. (64) we have to require that $c=u_{\tau}/\tau<1$, otherwise the
matching condition for the action cannot be met since the branch $\theta_{+}$
has a minimum at finite time (see Fig. 13).

As for the value of the action along the lower branch, Eq. (52) yields, if we
observe that $\sigma_{-}$ is of order $O(1)$,

$S_{-}=2(1+H-)\tau+2u_{\tau}+O(1)=\tau\left(\frac{u_{\tau}}{\tau}+1\right)^{2}+O(1),$
(65)

where we have used Eq. (63) as well. For the action along the upper branch
(50) we have to take into account that $\sigma_{+}$ can be expressed as
$-\ln(-2H_{+})+O(1)$ so that

$\displaystyle S_{+}$ $\displaystyle=$ $\displaystyle
2\tau+2u_{\tau}+2\ln(-2H_{+})+O(1)=4u_{\tau}+O(1),$ (66)

where we have employed Eq. (64). The kink’s position is determined by equating
the expressions (65) and (66), resulting in $u_{\tau}=\tau$, so that the kink
moves with unit speed. Finally, the minimized action to lowest order reads

$S=\left\\{\begin{array}[]{lcl}4u_{\tau}&\mbox{ for }u_{\tau}<\tau\\\
\tau(u_{\tau}/\tau+1)^{2}&\mbox{ for }u_{\tau}>\tau,\end{array}\right.$ (67)

which indeed describes the numerical findings shown in Fig. 7.

##  V Conclusion

We have studied in this paper a piecewise-constant SDE in order to understand
to what extent saddle-point approximations of the path integral representation
of the propagator, which is the basis of weak-noise approximations, can be
carried out. The advantage of the SDE that we have considered is that its
propagator is known exactly, so that any saddle-point approximations can be
compared to the exact result.

For this model, we have seen that saddle-point approximations are able to
reproduce some features of the propagator, such as its tail behavior and its
convergence towards the stationary distribution, but not the bimodality of the
exact propagator appearing at intermediate times. We have also seen that the
lowest-order approximation yields the correct large-deviation approximation of
the propagator, implying that the low-noise large-deviation theory of
Freidlin-Wentzell can be applied to this singular SDE. However, the
construction of higher-order corrections to this approximation is plagued by a
fundamental singularity of the Jacobian term of the path integral.

To remove this singularity, we have regularized the discontinuity of the SDE
with a smooth nonlinear drift involving a small parameter controlling the
limit to the piecewise smooth drift. The price paid for introducing this
regularization is that the two limits, small diffusion and small
regularization, do not commute and that the propagator of the regularized SDE
is no longer known exactly. Nevertheless, we have shown with simulation
results that the low-noise limit of the regularized SDE captures the main
features of the piecewise smooth SDE. In particular, the orbit or optimal path
structure of the regularized SDE in terms of direct, indirect and intermediate
paths is similar to the optimal path structure inferred heuristically for the
piecewise smooth SDE. In addition, the analysis of the regularized SDE
justifies the heuristic principles that we have defined and used to perform
the saddle-point approximation of the piecewise smooth SDE. For the
regularized SDE considered here, we have finally been able to study the
quasipotential associated with the propagator in a largely analytical way.
This is one of the few models for which such results can be obtained. In a
future study, we plan to consider higher-dimensional systems, such as models
in which inertia is present.

###### Acknowledgements.

Y.C. was supported by the Chinese Scholarship Council, the Hunan Provincial
Innovation Foundation for Postgraduates (Grant No. CX2011B011), and NUDT’s
Innovation Foundation (Grant No. B110205). W.J. gratefully acknowledges
support from EPSRC through Grant No. EP/H04812X/1 and No. SFB910 and the kind
hospitality by Eckehard Schöll and his group during a stay at TU Berlin.

##  Appendix A Weak-noise approximation of path integrals

Path integral formulations of the propagator of SDEs and their expansion for
weak-noise are well established in the literature. To set the notation and to
keep the presentation self-contained, we summarize here the essential
features, following the formulation presented in Chaichian and Demichev
(2001). Let us consider the one-dimensional Langevin equation

$\dot{v}=-f(v)+\sqrt{D}\xi(t),$ (68)

where $f(v)$ is a smooth force and the Gaussian white noise $\xi(t)$ obeys Eq.
(2). We can write down the conditional probability $p(v_{t},t|v_{0},0)$ by
using the path integral formula Chaichian and Demichev (2001),

$p(v_{t},t|v_{0},0)=\int_{(v_{0},0)}^{(v_{t},t)}\mathcal{D}[v]J[v]\mbox{e}^{-\frac{1}{4D}\int_{0}^{t}(\dot{v}+f(v))^{2}\mathrm{d}s}$
(69)

where $\int\mathcal{D}[v]$ denotes the Wiener measure and the Jacobian term

$J[v]=\exp\left(\frac{1}{2}\int_{0}^{t}f^{\prime}(v)\mbox{d}s\right)$ (70)

originates from the transformation $\xi(t)\rightarrow v(t)$. Putting the two
exponentials together, we thus express the kernel of the path integral in
terms of the action

$S[v]=\int_{0}^{t}L(v(s),\dot{v}(s))\mbox{d}s=\int_{0}^{t}\left((\dot{v}(s)+f(v(s)))^{2}-2Df^{\prime}(v(s))\right)\mbox{d}s.$
(71)

All trajectories contribute to the path integral (69), but for small $D$, the
largest contribution will come from the trajectory with smallest action. At
lowest order in $D$, this contribution is found by minimizing the action

$S^{(0)}[v]=\int_{0}^{t}\left(\dot{v}(s)+f(v(s))\right)^{2}\mbox{d}s,$ (72)

which does not take the contribution of Jacobian into account since it is
multiplied by $D$. The corresponding boundary value problem determined by the
Euler-Lagrange equation reads

$\ddot{v}^{(0)}(s)=f(v^{(0)}(s))f^{\prime}(v^{(0)}(s)),\quad
v^{(0)}(0)=v_{0},\quad v^{(0)}(t)=v_{t}.$ (73)

Given the path minimizing the action (72), the leading order approximation of
the propagator is thus given by

$p^{(0)}(v_{t},t|v_{0},0)=N_{1}\exp(-S^{(0)}[v^{(0)}]/(4D)),$ (74)

where the (time dependent) normalization $N_{1}$ can be computed a posteriori.

To improve this approximation, one can construct asymptotic series expansions
of the action, which do not, however, always ensure positivity of the
propagator. A simpler way to improve on the approximation (74) is to evaluate
the action of Eq. (71) with the Jacobian using the optimal path $v^{(0)}$; see
Eq. (73). This leads to the expression

$p^{(0;1)}(v_{t},t|v_{0},0)=N_{2}\exp(-S[v^{(0)}]/(4D)).$ (75)

A more coherent approach, which keeps the spirit of the saddle-point
approximation and ensures positivity of the propagator, is to evaluate the
minimizing path with the Jacobian, leading to the following Euler-Lagrange
equation:

$\ddot{v}^{(1)}(s)=f(v^{(1)}(s))f^{\prime}(v^{(1)}(s))-Df^{\prime\prime}(v^{(1)}(s)),\quad
v^{(1)}(0)=v_{0},\quad v^{(1)}(t)=v_{t}.$ (76)

In this case, the corresponding first order expression for the propagator is
given by

$p^{(1)}(v_{t},t|v_{0},0)=N_{3}\exp(-S[v^{(1)}]/(4D)).$ (77)

In general, it is not possible to solve the boundary value problem of the
Euler-Lagrange equations (73) or (76) analytically. Hence, numerical methods
such as the shooting method must be used. In this case, it is useful to
respect the underlying Hamiltonian structure of the Euler-Lagrange problem by
using symplectic integration methods which preserve the Hamiltonian,

$H^{(0)}(v,p_{v})=p_{v}^{2}/2-f^{2}(v)/2$ (78)

without the Jacobian or

$H(v,p_{v})=p_{v}^{2}/2-(f^{2}(v)+Df^{\prime}(v))/2$ (79)

with the Jacobian. In this case, a symplectic Euler scheme or the Störmer-
Verlet scheme, for example, can be applied to integrate the corresponding
canonical equations of motion.

##  Appendix B Energy-time relation

###  B.1 Monotonicity of $\theta_{-}$

Consider the expression defined in Eq. (42) either for $u_{\tau}<0<u_{0}$ and
$H>0$ or for $0<u_{\tau}<u_{0}$ and $H>H_{min}=-\tanh^{2}(u_{\tau})/2$. The
argument of the logarithm in Eq. (42) can be written as

$g(x)=\frac{a+\sqrt{x+a^{2}}}{b+\sqrt{x+b^{2}}}$ (80)

and it is easy to see that $g(x)>1$ in the given parameter ranges, i.e, either
$b<0<a$ and $x>0$, or $0<b<a$ and $x>-b^{2}$. In addition, it follows by
differentiation that $g^{\prime}(x)<0$, i.e., $g$ is monotonically decreasing.
Hence $\theta_{-}$ considered as a function of $H$ is the product of two
positive monotonically decreasing functions so that $\theta_{-}$ itself is
monotonically decreasing.

###  B.2 Critical points of $\theta_{+}$

Since Eq. (47) is a symmetric expression in the first two arguments it is
sufficient to consider the case $0<u_{\tau}<u_{0}$ and
$H_{min}=-\tanh^{2}(u_{\tau})/2<H<0$. The critical points are determined by
the vanishing derivative of $\theta_{+}$. It is more convenient to consider
the expression in terms of the new variable $\chi\in[0,u_{\tau}]$ defined by
$H=-\tanh^{2}(\chi)/2$, where $-\chi$ represents the turning point of the
indirect path, introduced in Sec. IV.1.2. Then differentiation gives

$\displaystyle\frac{\partial\theta_{+}(u_{0},u_{\tau},-\tanh^{2}(\chi)/2)}{\partial\chi}=h(\sinh(u_{0}),\sinh(\chi))+h(\sinh(u_{\tau}),\sinh(\chi)),$
(81)

where we have introduced

$h(a,z)=z\ln\left(a/z+\sqrt{(a/z)^{2}-1}\right)-\frac{(1+z^{2})a}{z\sqrt{a^{2}-z^{2}}},\quad(0<z<a).$
(82)

The critical points of $\theta_{+}$ are thus determined by the solutions of
the equation

$h(\sinh(u_{0}),z)+h(\sinh(u_{\tau}),z)=0.$ (83)

We now show that the right-hand side of Eq. (81) is a convex function of
$z=\sinh(\chi)$, so there exist at most two solutions. To do so, compute the
second derivative of Eq. (82):

$\frac{\partial h(a,z)}{\partial
z}=\ln\left(a/z+\sqrt{(a/z)^{2}-1}\right)+a\frac{z^{4}-2(a^{2}+1)z^{2}+a^{2}}{z^{2}\left(a^{2}-z^{2}\right)^{3/2}}$
(84)

and

$\frac{\partial^{2}h(a,z)}{\partial
z^{2}}=-a\frac{z^{4}(6+2a^{2})+z^{2}(-5a^{2}+a^{4})+2a^{4}}{z^{3}\left(a^{2}-z^{2}\right)^{5/2}}.$
(85)

It is easy to see that the numerator in Eq. (85) is positive by completing the
square $6(z^{2}-a^{2}/2)^{2}$. Hence, the second derivative is negative for
any (positive) values of $a$. The same holds for Eq. (83), proving convexity.

Finally, it is obvious from Eq. (47) that $\theta_{+}(u_{0},u_{\tau},H)$ tends
to $\infty$ as $H\uparrow 0$ and that $\theta_{+}(u_{0},u_{\tau},H_{min})$ is
finite and positive, resulting altogether in a cubic shape for the graph of
$\theta_{+}(u_{0},u_{\tau},H)$.

## References

  * Filippov (1988) A. F. Filippov, _Differential Equations with Discontinuous Righthand Sides_ (Kluwer, Dordrecht, 1988).
  * di Bernardo et al. (2008) M. di Bernardo, C. Budd, A. Champneys, and P. Kowalczyk, _Piecewise-Smooth Dynamical Systems: Theory and Applications_ (Springer, Berlin, 2008).
  * Makarenkov and Lamb (2012) O. Makarenkov and J. S. W. Lamb, Physica D 241, 1826 (2012).
  * Biemond et al. (2012) B. J. J. Biemond, N. van de Wouw, and H. Nijmeijer, Physica D 241, 1882 (2012).
  * Colombo et al. (2012) A. Colombo, M. di Bernardo, S. J. Hogan, and M. R. Jeffrey, Physica D 241, 1845 (2012).
  * Bowden and Tabor (1950) F. P. Bowden and D. Tabor, _Friction and Lubrication of Solids_ (Oxford University Press, Oxford, 1950).
  * Elmer (1997) F.-J. Elmer, J. Phys. A: Math. Gen. 30, 6057 (1997).
  * Persson (1998) B. N. J. Persson, _Sliding Friction: Physical Principles and Applications_ (Springer, Berlin, 1998).
  * Berger (2002) E. J. Berger, Appl. Mech. Rev. 55, 535 (2002).
  * Jeffrey (2011) M. R. Jeffrey, Phys. Rev. Lett. 106, 254103 (2011).
  * Reimann (2002) P. Reimann, Phys. Rep. 361, 57 (2002).
  * Caughey and Dienes (1961) T. K. Caughey and J. K. Dienes, J. Appl. Phys. 32, 2476 (1961).
  * Atkinson and Caughey (1968) J. D. Atkinson and T. K. Caughey, Int. J. Nonlin. Mech. 3, 399 (1968).
  * de Gennes (2005) P.-G. de Gennes, J. Stat. Phys. 119, 953 (2005).
  * Hayakawa (2005) H. Hayakawa, Physica D 205, 48 (2005).
  * Baule and Sollich (2012) A. Baule and P. Sollich, Europhys. Lett. 97, 20001 (2012).
  * Touchette et al. (2010) H. Touchette, E. V. der Straeten, and W. Just, J. Phys. A: Math. Theor. 43, 445002 (2010).
  * Chaudhury and Mettu (2008) M. K. Chaudhury and S. Mettu, Langmuir 24, 6128 (2008).
  * Goohpattader et al. (2009) P. S. Goohpattader, S. Mettu, and M. K. Chaudhury, Langmuir 25, 9969 (2009).
  * Goohpattader and Chaudhury (2010) P. S. Goohpattader and M. K. Chaudhury, J. Chem. Phys. 133, 024702 (2010).
  * Gnoli et al. (2013) A. Gnoli, A. Puglisi, and H. Touchette, Europhys. Lett. 102, 14002 (2013).
  * Simpson and Kuske (2012) D. J. W. Simpson and R. Kuske, arXiv:1204.5792 (2012).
  * Karatzas and Shreve (1984) I. Karatzas and S. E. Shreve, Ann. Prob. 12, 819 (1984).
  * Zhang (1990) W. Zhang, IEEE Trans. Automat. Contr. 35, 980 (1990).
  * Benes et al. (1980) V. E. Benes, L. A. Shepp, and H. S. Witsenhausen, Stochastics 4, 39 (1980).
  * Browne and Whitt (1995) S. Browne and W. Whitt, in _Advances in Queueing_ , edited by J. Dshalalow (CRC Press, Boca Raton, FL, 1995), pp. 463–480.
  * Gradinaru et al. (2001) M. Gradinaru, S. Herrmann, and B. Roynette, Ann. Inst. Poincaré B 37, 555 (2001).
  * Korostelev and Leonov (1993) A. P. Korostelev and S. L. Leonov, Prob. Theory Relat. Fields 94, 317 (1993).
  * Chiang and Sheu (2000) T.-S. Chiang and S.-J. Sheu, Ann. Prob. 28, 140 (2000).
  * Boué et al. (2000) M. Boué, P. Dupuis, and R. S. Ellis, Prob. Theory Relat. Fields 116, 125 (2000).
  * Chaichian and Demichev (2001) M. Chaichian and A. Demichev, _Path Pntegrals in Physics: Stochastic Processes and Quantum Mechanics_ (IOP, London, 2001), Vol. 1.
  * Graham (1989) R. Graham, in _Noise in Nonlinear Dynamical Systems_ , edited by F. Moss and P. V. E. McClintock (Cambridge University Press, Cambridge, 1989), Vol. 1, pp. 225–278.
  * Touchette (2009) H. Touchette, Phys. Rep. 478, 1 (2009).
  * Baule et al. (2010) A. Baule, E. G. D. Cohen, and H. Touchette, J. Phys. A: Math. Theor. 43, 025003 (2010).
  * Baule et al. (2011) A. Baule, H. Touchette, and E. G. D. Cohen, Nonlinearity 24, 351 (2011).
  * Freidlin and Wentzell (1984) M. I. Freidlin and A. D. Wentzell, _Random Perturbations of Dynamical Systems_ , Grundlehren der Mathematischen Wissenschaften Vol. 260 (Springer, New York, 1984).
  * Kloeden and Platen (1992) P. Kloeden and E. Platen, _Numerical Solution of Stochastic Differential Equations_ (Springer, Berlin, 1992).
  * Graham and Tél (1985) R. Graham and T. Tél, Phys. Rev. A 31, 1109 (1985).
  * Cassing and Knoll (1981) W. Cassing and J. Knoll, Phys. Lett. B 106, 443 (1981).