Sticky Brownian Motion and its Numerical Solution

TL;DR

This paper introduces a new, efficient simulation method for sticky Brownian motion, a diffusion process with boundary interaction, with applications in biology, materials science, and finance, and demonstrates significant speed improvements over existing methods.

Contribution

It presents a simple sticky random walk model for simulating sticky Brownian motion, offering substantial computational speed gains and insights into its properties.

Findings

Sticky random walk is 100 to 10,000 times faster than existing methods.

The model accurately captures mesoscale particle behavior near boundaries.

Potential extension to multi-dimensional sticky diffusions is discussed.

Abstract

Sticky Brownian motion is the simplest example of a diffusion process that can spend finite time both in the interior of a domain and on its boundary. It arises in various applications such as in biology, materials science, and finance. This article spotlights the unusual behavior of sticky Brownian motions from the perspective of applied mathematics, and provides tools to efficiently simulate them. We show that a sticky Brownian motion arises naturally for a particle diffusing on $\mathbb{R}_+$ with a strong, short-ranged potential energy near the origin. This is a limit that accurately models mesoscale particles, those with diameters $\approx 100$nm-$10\mu$m, which form the building blocks for many common materials. We introduce a simple and intuitive sticky random walk to simulate sticky Brownian motion, that also gives insight into its unusual properties. In parameter regimes of practical interest, we show this sticky random walk is two to five orders of magnitude faster than alternative methods to simulate a sticky Brownian motion. We outline possible steps to extend this method towards simulating multi-dimensional sticky diffusions.

Equations

$$d{X}_t^{ϵ}=-{\partial }_x{U}^{ϵ}(X_t^{ϵ})dt+\sqrt{2}d{W}_t,$$

$$\underset{ϵ\to 0}{\lim }{\int }_0^{ϵ}{e}^{-U^{ϵ}(x)}dx=\kappa .$$

$$\sqrt{2\pi }\underset{ϵ\to 0}{\lim }\frac{ϵe^{-U^{ϵ}(x_0)}}{\sqrt{{\partial }_{xx}{U}^{ϵ}(x_0)}}=\kappa ,$$

$$U^{ϵ}(x)=D_e(1-e^{-a(x-x_0)}{)}^2-D_e$$

$$U^{ϵ}(x)=D_e\left({\left(\frac{1}{x-x_0+1}\right)}^{2m}-2{\left(\frac{1}{x-x_0+1}\right)}^m\right)$$

$$\underset{ϵ\to 0}{\lim }{\tau }^{ϵ}(0)=\kappa \ell +\frac{{\ell }^2}{2}.$$

$${\partial }_x(e^{-U^{ϵ}}{\partial }_x{\tau }^{ϵ})=-e^{-U^{ϵ}}\text{on [0,ℓ] with b.c. ∂xτϵ(0)=0 and τϵ(ℓ)=0.}$$

$${\tau }^{ϵ}(x)={\int }_x^{\ell }{\int }_0^r{e}^{U^{ϵ}(r)-U^{ϵ}(s)}dsdr$$

$${\tau }^{ϵ}(0)$$

$$+{\int }_{ϵ}^{\ell }{e}^{-U^{ϵ}(s)}{\int }_s^{\ell }{e}^{U^{ϵ}(r)}drds.$$

$$\left|{\int }_0^{ϵ}{e}^{-U^{ϵ}(s)}{\int }_s^{ϵ}{e}^{U^{ϵ}(r)}drds\right|\le {\int }_0^{ϵ}{\int }_s^{ϵ}{e}^{0\vee (-U(s)+u_{\ast })}drds\to 0\text{ as }ϵ\to 0.$$

$${\partial }_t{p}^{ϵ}(x,t)+{\partial }_x{j}^{ϵ}(x,t)=0,$$

$${\tilde{X}}_{k+1}=|{\tilde{X}}_k-{\partial }_x{U}^{ϵ}({\tilde{X}}_k)\delta t+\sqrt{2}(W(t_{k+1})-W(t_k))|,{\tilde{X}}_0=x\ge 0.$$

$$\rho (x,t)=p(x,t)(1+\kappa \delta (x)),$$

$${\partial }_{t}p={\partial }_{xx}p,\text{with b.c. }\kappa {\partial }_{xx}p={\partial }_{x}p\text{ at }x=0,$$

$$p^{ϵ}(x,t)\sim \{\begin{array}{ll}{p}_0(x,t)+p_1(x,t)+\cdots ,& x>ϵ,\\ q_0(x,t)+q_1(x,t)+\cdots ,& x\le ϵ,\end{array}$$

$${\partial }_t{p}_0={\partial }_{xx}{p}_0,x\ge ϵ,$$

$$0={\partial }_X({\partial }_X{U}^{ϵ}(ϵX)q_0)+{\partial }_{XX}{q}_0,X\le 1$$

$$q_0(ϵX,t)=a(t)e^{-U^{ϵ}(ϵX)},X\le O(1)$$

$$p_0(ϵ,t)=q_0(ϵ,t)⟹a(t)=p_0(0,t)$$

$$\frac{d}{dt}\left({\int }_0^{ϵ}{q}_0(x,t)dx+{\int }_{ϵ}^{\infty }{p}_0(x,t)dx\right)=0⟹\kappa a^{\prime }(t)-{\partial }_x{p}_0(0,t)=0,$$

$$\kappa {\partial }_t{p}_0{|}_{x=0}={\partial }_x{p}_0{|}_{x=0}⟺\kappa {\partial }_{xx}{p}_0{|}_{x=0}={\partial }_x{p}_0{|}_{x=0}.$$

$$\pi (x)=Z^{-1}(1+\kappa \delta (x)),\text{where}Z=\kappa +L.$$

$$\rho (x,t)={\rho }_0(t){\mu }_0+{\rho }_1(x,t){\mu }_1,$$

$${\rho }_0(t)=\kappa {\rho }_1(0,t)$$

$$\begin{array}[]{rll}\partial_{t}\rho_{0}=&\partial_{x}\rho_{1}&\text{for }x=0\,,\\ \partial_{t}\rho_{1}=&\partial_{xx}\rho_{1}&\text{for }x\in(0,\infty)\,.\end{array}$$

$$⟨f,{\partial }_t\rho ⟩=⟨\mathcal{L}f,\rho ⟩.$$

$$⟨f,{\partial }_t\rho ⟩$$

$$=(f\kappa {\partial }_{xx}p){|}_{x=0}-(f{\partial }_{x}p){|}_{x=0}+(p{\partial }_{x}f){|}_{x=0}+{\int }_0^{\infty }p{\partial }_{xx}fdx$$

$$=(p{\partial }_{x}f){|}_{x=0}-(\kappa p{\partial }_{xx}f){|}_{x=0}+\int \rho {\partial }_{xx}fdx.$$

$$\mathcal{L}f={\partial }_{xx}f,\text{with b.c. }\kappa {\partial }_{xx}f={\partial }_{x}f\text{ at }x=0.$$

$${\partial }_{xx}\tau (x)=-1\text{on [0,ℓ] with b.c. κ∂xxτ(0)=∂xτ(0) and τ(ℓ)=0}.$$

Full text

{AMS}

60H10, 65C30 (60J60, 60J65, 35K05, 35K20, 65M06)

1 Introduction

Sticky diffusion processes are solutions to stochastic differential equations (SDE) which can ‘stick’ to, i.e. spend finite time on, a lower-dimensional boundary. The sticking is reversible, so the process can hit the boundary and leave again, and while on the boundary it can move according to dynamics that are different from those in the interior, even when continuously extended to the boundary. A simple example is a (root-2) Brownian motion which can stick to the origin, called a sticky Brownian motion, whose forward and backward Kolmogorov equations are identically $\partial_{t}f=\partial_{xx}f$ with boundary condition $\partial_{x}f|_{x=0}=\kappa\partial_{xx}f|_{x=0}$, where $\kappa\geq 0$ is a parameter measuring how sticky the boundary is. On the other hand Dirichlet ($f|_{x=0}=0$), Neumann ($\partial_{x}f|_{x=0}=0$), Wentzell ($\partial_{xx}f|_{x=0}=0$) and Robin ($\partial_{x}f|_{x=0}=\kappa f_{x=0}$) boundary conditions lead to stopped, reflected, absorbed and elastic Brownian motions, respectively [48].

Discovered in the 1950s by Feller in an attempt to find the most general behaviour of one-dimensional diffusion processes at a boundary [27, 80], sticky diffusions have been studied in the theoretical mathematical literature for several decades. Probabilists have studied the construction and properties of sticky diffusions via martingales or other representations such as random walks in random environments or interacting particle systems [28, 29, 87, 45, 86, 53, 44, 1, 91, 75, 76, 5], and very recently, sticky diffusions have been used to expand the scope of probabilistic coupling techniques [43, 20, 96]. In PDE theory, analysts have studied well-posedness and semigroup generation for parabolic and elliptic problems with sticky boundary conditions, called generalized Wentzell or Wentzell-Robin boundary conditions in the PDE literature [67, 93, 2, 3, 25, 24].

In applications, sticky diffusions arise in a variety of models of physical and natural processes which are naturally described by a set of variables that can change dimension. Examples arise in biology, where molecules diffuse near a sticky wall or cell membrane [31, 37]; in epidemics, [10], where the concentration of pathogens in an individual can be sticky at concentration level zero; in operations research, as a particular limit of storage processes modeling queues, inventories, insurance risks, etc [38]; and in finance, where the evolution of interest rates can be sticky near zero [65, 46].

Our own interest is in the dynamics of mesoscale particles, those with diameters of $\approx 100nm-10\mu m$, which occur widely and form the building blocks of common materials like paint, toothpaste, concrete, ketchup, and many others [66]. Such particles interact attractively over ranges typically much smaller than their diameters [70, 41]. Remarkably, systems with short-ranged interactions are often insensitive to the exact shape of the attractive well of the interaction potential, with most behaviour depending only on one or two parameters such as some combination of the well depth and well width [77, 81]. Therefore, it is often effective to model such systems in the sticky limit, where the well width is taken to zero, and the well depth to infinity, such that the probability of forming a contact remains constant. In this limit the dynamics of the collection of particles approaches a sticky diffusion process, with a boundary when a pair of particles are exactly in contact [6, 85, 23, 72, 42, 47]. The sticky boundary conditions behave in a similar way to holonomic constraints in molecular dynamics, which eliminate fast bond-length or bond-angle vibrations, and thus reveal the molecular structure more clearly and allow a larger time step in simulations [84].

In the sticky limit, one may be interested in studying the forward and backward Kolomgorov equations to obtain analytical, asymptotic, or numerical solutions that give physical insight, or in simulating the sticky processes themselves to obtain pathwise results. However, neither of these goals is readily attainable: the first, since sticky processes are relatively unknown in the applied math community, techniques to study them are rare and usually invented on a case-by-case basis. Indeed, by-and-large, applied mathematicians focus on PDEs with classical boundary conditions, and sticky diffusions are beyond this scope, since their transition probability measures have a part that is singular with respect to the Lebesgue measure in the domain. The second, because there are currently no methods to simulate a sticky diffusion directly: there is no practical way to extend existing methods for discretizing SDEs based on choosing discrete time steps, such as Euler-Maruyama or its variants [50, 35, 7], to sticky processes; a rough explanation for why is that in these methods one will never hit the boundary exactly. One can approximate a sticky diffusion by a reflecting diffusion with an artificial force near the boundary to encourage the particle to stay there when it gets near, but for a good approximation, the force must be strong and short-ranged. Most SDE solvers are explicit, especially in molecular dynamics applications where evaluating forces is the most costly step, so one must take a prohibitively small timestep to resolve these forces, which unfortunately, severely limits the timescales one can simulate.

Our aim in this article is twofold: one, we wish to bring the topic of sticky diffusions and their associated PDEs to the attention of the applied math community, and to introduce tools to study them from an applied math perspective. To this end, we show how a sticky Brownian motion arises as a limit of reflected Brownian motions with a strong short-ranged force at the origin, and discuss the limiting forward and backward equations, which must be used with care because of their unusual boundary conditions (Section 2.) Two, we wish to introduce a numerical method to simulate a sticky diffusion, which simulates the process directly without introducing an artificial force, and which allows one to take a relatively large timestep. The method is based on discretizing the increments of the process in space, rather than in time, and constructing a Markov jump process whose generator locally approximates the generator of the sticky diffusion (Section 3.) We derive a Feynman-Kac formula to show that this Markov jump process can be used to solve certain PDEs to second-order accuracy in the spatial step. Basic implementations of the two main numerical approximations used in the paper are provided in Appendix A.

Throughout the paper we focus on a one-dimensional sticky Brownian motion, because this illustrates most of the key ideas and differences from traditional diffusion processes with a minimum of technical difficulties. We discuss the connection between our approaches to understanding sticky Brownian motion and those taken in the earlier probability literature (Section 4.) We expect the methods we introduce to be fully adaptable to higher-dimensional diffusion processes, and outline the steps required to do so in the conclusion (Section 5.)

2 Sticky Brownian motion

In this section we consider how a sticky Brownian motion (SBM) arises naturally for a particle diffusing in a potential energy landscape that has a strong, short-ranged potential well near the origin.111In Sections 2-3.4, what we call a sticky Brownian motion is traditionally called a root-2 sticky Brownian motion, because it is scaled by a factor of $\sqrt{2}$ from a traditional Brownian motion. In Section 4 we use SBM to refer to a traditional sticky Brownian motion. In turn, this motivates the generator and the backward and forward Kolmogorov equations for an SBM, equations and operators which require working in an unusual function space because the transition probabilities will have a singular part.

2.1 Setup

Consider a diffusion process $X^{\epsilon}_{t}$ on $\mathbb{R}_{\geq 0}$, depending on a parameter $\epsilon>0$, which solves

[TABLE]

with a reflecting boundary condition at the origin. Here $W_{t}$ is a standard Brownian motion, and $U^{\epsilon}(x):\mathbb{R}\to\mathbb{R}$ is a function parameterized by $\epsilon$. If $X^{\epsilon}_{t}$ is the position at time $t$ of a particle moving on the real axis, then the force it feels is $-\partial_{x}U^{\epsilon}(X^{\epsilon}_{t})$. Recall that, if $\int_{0}^{\infty}e^{-U^{\epsilon}(x)}dx<\infty$, then $e^{-U^{\epsilon}(x)}$ is the non-normalized stationary (or equilibrium) probability density of $X^{\epsilon}_{t}$.

Equation (1) is a special case of the Brownian dynamics equations, which are a good model for the dynamics of mesoscale particles in a fluid [30]. For a more general system of particles, $X_{t}^{\epsilon}$ would be the configuration, a vector of particle positions, $U^{\epsilon}(x)$ would represent the potential energy of a particular configuration (nondimensionalized by temperature), which is usually a sum of the potential energy between each pair of particles, and the equations may additionally include a friction tensor, depending on the configuration and the velocity, modeling hydrodynamic interactions between particles. We specialize to a single scalar equation and ignore friction, but still think of $U^{\epsilon}(x)$ as a potential energy, representing the energy of the particle $X_{t}^{\epsilon}$ as a function of its distance to another particle. Physically, this equation could be realized by holding one particle in place at the origin while another is particle moves on a line.

An important point is that for mesoscale particles, $U^{\epsilon}(x)\approx 0$ outside an interval that is very narrow compared to the particles’ diameters [70, 41]. For example, for particles interacting with a so-called “depletion” interaction [4], the range of the interaction in a typical experiment was estimated to be about 5% of the particles’ diameters [71]. Therefore, when two particles’ surfaces are within this range of each other, they stay close together for a long time, but when they are further apart they don’t feel each other at all, and diffuse independently. An even shorter range is achieved by particles that interact via sticky single-stranded DNA coated on their surfaces: here, the range of the interaction is the average radius of a coiled DNA strand, which is typically around 10nm, about 1% of the diameter of a 1$\mu$m particle [90]. Such DNA-coated particles are studied extensively because the DNA allows one to code complex interactions between different types of particles, and hence, to program them to assemble into a great many different materials [68, 12, 90, 83, 95].

2.2 Assumptions on the potential energy function

With these remarks in mind, we consider a family of potential energy functions $(U^{\epsilon}(x))_{\epsilon>0}$ with a narrow, deep attractive well, which becomes narrower and deeper as $\epsilon\to 0$. We consider the dynamics of $X^{\epsilon}_{t}$ as $\epsilon\to 0$, and call this the sticky limit. Specifically, we impose the following assumptions on $(U^{\epsilon}(x))_{\epsilon>0}$.

{assumption}

For any $\epsilon>0$, $U^{\epsilon}$ is a function in $C^{2}(\mathbb{R})$ satisfying:

• (A1)

$U^{\epsilon}(x)$, $\partial_{x}U^{\epsilon}(x)$, $\partial_{xx}U^{\epsilon}(x)\leq O(\epsilon)$ for $x\geq\epsilon$;

• (A2)

$U^{\epsilon}(x)$ possesses a unique local minimum in $(0,c\epsilon)$ with minimizer $x_{0}$ and no local maximum for $x>0$;

• (A3)

There exists $\kappa\geq 0$ such that

[TABLE]

We briefly comment on the physical interpretation of (A1)-(A3). Assumption (A1) ensures the potential and its first two derivatives are negligible outside of the interval $(0,\epsilon)$, which we call the boundary layer. Outside the boundary layer, $X^{\epsilon}_{t}$ feels virtually no force and simply diffuses.

Assumption (A2) ensures that the dynamics in (1) has at most one timescale within the boundary layer. A typical timescale for a diffusion process is its mean first passage time (MFPT) to overcome an energy barrier $\Delta U$. Under Assumption (A2), the barrier to leaving the interval $(0,\epsilon)$ is the depth of the potential at its minimum, leading to an MFPT of approximately $\epsilon e^{\Delta U}$ [32]. The assumption rules out pathological potentials with many large oscillations in the boundary layer, which would give rise to longer dynamical timescales in the boundary layer.

Assumption (A3) is the one that gives rise to stickiness at the origin. It requires that the measure of $(0,\epsilon)$ with respect to the weighted Lebesgue measure $e^{-U(x)}dx$, or in the language of physics, the partition function for this interval, approaches a constant. This constant $\kappa$ determines how sticky the origin is – larger $\kappa$ means the process will spend longer near the origin on average. For this reason we call $\kappa$ the sticky parameter. Applying Laplace asymptotics to (2) shows that

[TABLE]

where $x_{0}$ is the minimizer of $U^{\epsilon}(x)$ in $(0,\epsilon)$. Ignoring the second derivative shows that the depth must scale very nearly as $|U^{\epsilon}(x_{0})|\sim|\log\epsilon|$, the logarithm of the width of the potential. This implies the timescale computed from the MFPT for leaving the boundary layer is $\epsilon e^{\Delta U}\sim O(1)$. If the scaling of the depth of the well is larger than $|\log\epsilon|$, then large deviation theory would be more appropriate to describe the dynamics of the limiting process [14, 36]. However, if the scaling is smaller, then the limiting process spends no time on the boundary, and simply reflects off of it, as we will see momentarily.

The requirement that $U^{\epsilon}$ has two derivatives is not necessary, but is included to simplify our calculations and avoid dealing with discontinuities in the coefficients of (1). An example of a potential energy function which doesn’t satisfy this condition, but which is commonly used to model short-ranged potentials, is the square-well potential.222Square-well potentials have the form $U^{\epsilon}(x)=C_{\epsilon}$ for $x\in[0,\epsilon]$, $U^{\epsilon}(x)=0$ for $x>\epsilon$, $U^{\epsilon}(x)=\infty$ for $x<0$; where $C_{\epsilon}$ is a constant which depends on $\epsilon$. By introducing a smooth approximation, we expect our asymptotic results to hold for a square-well potential as well, though we do not pursue this here.

Assumption 2.2 can be verified for two families of potentials frequently used to model attractive interactions between mesoscale particles, the Morse and generalized Lennard-Jones potential energy functions [18, 69, 72, 88, 11, 94].

Example 2.1** (Morse Potential).**

Fix $\kappa\geq 0$, let $\epsilon>0$, and consider the potential energy function defined by

[TABLE]

*with parameters $a=1/\epsilon^{2}$, $x_{0}=\epsilon^{2}$, and $D_{e}$ defined implicitly via $e^{D_{e}}\frac{\sqrt{\pi}}{a\sqrt{D_{e}}}=\kappa$. This potential is illustrated in Figure 1 for several different values of $\epsilon$. It has a unique global minimum at $x_{0}\in(0,\epsilon)$ with depth $U^{\epsilon}(x_{0})=-D_{e}$, and hence, (A2) holds. The width of the basin of attraction of the minimum is $O(a^{-1})=O(\epsilon^{2})$. Moreover $U^{\epsilon}(x)$ and all of its derivatives are exponentially small for $x\geq\epsilon$, which implies that (A1) holds. Lastly, as $\epsilon\to 0$, $D_{e}\to\infty$ with parameters chosen as above, straightforward Laplace asymptotics gives that $\lim_{\epsilon\to 0}\dfrac{\int_{0}^{\epsilon}e^{-U^{\epsilon}(x)}dx}{e^{-U^{\epsilon}(x_{0})}\frac{\sqrt{2\pi}}{\sqrt{(U^{\epsilon})^{\prime\prime}(x_{0})}}}=\lim_{\epsilon\to 0}\dfrac{\int_{0}^{\epsilon}e^{-U^{\epsilon}(x)}dx}{e^{D_{e}}\frac{\sqrt{\pi}}{a\sqrt{D_{e}}}}=1\,.$ Thus, (A3) holds. *

Example 2.2** (Lennard-Jones($2m$,$m$) potential).**

Fix $\kappa\geq 0$, let $\epsilon>0$, and consider the potential energy function defined by

[TABLE]

*with parameters $m=1/\epsilon^{2}$, $x_{0}=\epsilon^{2}$ and $D_{e}$ defined implicitly via $e^{D_{e}}\frac{\sqrt{\pi}}{m\sqrt{D_{e}}}=\kappa$. This potential has a unique global minimum at $x_{0}=\epsilon^{2}$ with depth $U^{\epsilon}(x_{0})=-D_{e}$, and hence, (A2) holds. The width of the basin of attraction of the minimum is $O(m^{-1})=O(\epsilon^{2})$. Moreover $U^{\epsilon}(x)$ and all of its derivatives are exponentially small for $x\geq\epsilon$, which implies that (A1) holds. Similar to the preceding example, one can verify that (2) holds. Hence, Assumption 2.2 is satisfied. *

2.3 Dynamics of $X^{\epsilon}_{t}$ for $\epsilon\ll 1$

Consider the dynamics of $X^{\epsilon}_{t}$ in (1) when $\epsilon$ is small. When $X^{\epsilon}_{t}$ is far enough away from the origin, it feels no force, and simply diffuses, like a Brownian motion. When $X^{\epsilon}_{t}$ is within a distance of $\epsilon$ from the origin, it feels a strong force keeping it near the minimum of $U^{\epsilon}$ for a long time, until an occasional large fluctuation pushes it out of the range of the force.

How long does it stay near the origin, and does this time remain significant as $\epsilon\to 0$? We start by computing the mean first-passage time (MPFT) to escape from a region near the origin.

Lemma 2.3**.**

Let $\tau^{\epsilon}(x)=\mathbb{E}\left(\inf\{t\geq 0\mid X^{\epsilon}_{t}>\ell,~{}X_{0}=x\in[0,\ell]\}\right)$ be the MFPT of $X_{t}^{\epsilon}$ out of $[0,\ell]$ for some $\ell>0$ with initial condition $x\in[0,\ell]$. Then

[TABLE]

When $\kappa=0$, we recover the MFPT of a reflecting Brownian motion starting at [math]. The time scales as the distance squared, $\tau^{\epsilon}(0)\sim O(\ell^{2})$, a traditional diffusive scaling. When $\kappa\to\infty$, the MFPT is infinite, consistent with the MFPT of an absorbing Brownian motion starting at [math]. For intermediate $\kappa$, the MFPT at the origin scales as $\tau^{\epsilon}(0)\sim O(\ell)$, which is a ballistic scaling – slower (for small $\ell$) than the diffusive scaling. Therefore, we expect the limiting probability density near the origin to be correspondingly large.

Proof 2.4**.**

The MFPT $\tau^{\epsilon}(x)$ satisfies the boundary value problem

[TABLE]

By integrating twice, the semi-analytic solution to this equation is given by

[TABLE]

and at $x=0$ we obtain,

[TABLE]

*Assumptions 2.2 (A1) and (A3) imply that the last two terms converge to the right-hand side of (3), while Assumption 2.2 (A3) implies the first term converges to 0.333 To see this last point, let $x_{0}$ be the minimizer of $U(s)$ on $[0,\epsilon]$, and let $u_{*}=\max_{s\in[x_{0},\epsilon]}U(s)$ be the maximum of the energy to the right of the minimum. By Assumption 2.2 (A2), there is at most one point $x_{*}\in[0,x_{0})$ such that $U(x_{*})=u_{*}$, and since the energy monotonically increases as $s$ decreases below $x_{*}$, we have $-(U(s)-U(r))\leq 0$ for $r\geq s$, $s\leq x_{*}$. For $s>x_{*}$, $r\geq s$, we have $-(U(s)-U(r))\leq-U(s)+u_{*}$. Therefore by Assumption 2.2 (A3),

$\left|\int_{0}^{\epsilon}e^{-U^{\epsilon}(s)}\int_{s}^{\epsilon}e^{U^{\epsilon}(r)}drds\right|\leq\int_{0}^{\epsilon}\int_{s}^{\epsilon}e^{0\vee(-U(s)+u_{*})}drds\to 0\quad\text{ as }\epsilon\to 0\,.$

To gain further insight into the dynamics, consider the evolution of the probability density $p^{\epsilon}(x,t)$ of $X^{\epsilon}_{t}$. It evolves according to the Fokker-Planck equation

[TABLE]

where $j^{\epsilon}(x,t)=-\partial_{x}U^{\epsilon}(x)p^{\epsilon}(x,t)-\partial_{x}p^{\epsilon}(x,t)$ is the associated probability current, or probability flux. The equation is to be solved in the domain $x>0$, $t>0$ with boundary condition $j^{\epsilon}(0,t)=0$ for all $t\geq 0$, and given initial condition $p^{\epsilon}(x,0)=\phi(x)$.

For this one-dimensional problem we could compute $p^{\epsilon}(x,t)$ to high accuracy by solving the PDE (4) numerically. However, such an approach won’t work in the high dimensions characteristic of systems of particles, so to illustrate the difficulties that may arise, we simulate $10^{6}$ trajectories of (1) numerically and estimate the density by kernel density estimation. We use a Morse potential for $U^{\epsilon}$ as in Example 2.1, and a symmetrized Euler-Maruyama approximation to solve (1), a method that is first-order weakly accurate method for domains with a smooth enough boundary [35, 7]. The method works as follows: given a time-step size $\delta t$, we set $t_{k}=k\delta t$ for $k\in\mathbb{N}_{\geq 0}$, and compute the approximation $\tilde{X}_{k}$ to $X^{\epsilon}(t_{k})$ as

[TABLE]

See Appendix A for Matlab code implementing this method. Unfortunately, when the Morse potential used is short-ranged and strong, an accurate approximation requires a small time step size: when $\kappa=1,a\approx 118$ the time step size required for a qualitatively correct solution is approximately $\delta t\leq 2.5\times 10^{-6}$ – Figure 2 shows that larger timesteps fail to see the boundary layer entirely.

Figure 3 plots $p^{\epsilon}(x,1)$ on a log scale for various values of $\epsilon$, and shows that as $\epsilon$ decreases, $p^{\epsilon}(x,1)$ has an increasingly large peak in density near the origin. Although the width of the peak decreases, the total probability in the peak remains nearly constant with $\epsilon$. Therefore, we expect the probability density to contain a singularity of the form $\delta(x)$ when $\epsilon\to 0$.

Furthermore, notice that $p^{\epsilon}(x,1)$ is slowly varying for $x\geq\epsilon$ and rapidly varying for $0\leq x\leq\epsilon$. This suggests using the method of matched asymptotic expansions to study (4) as $\epsilon\to 0$ [40, 49]. This approach, which is very similar to boundary-layer theory in fluids, proceeds by finding a local ‘inner’ solution near the origin and a local ‘outer’ solution far enough away from the origin, and then matching these local solutions in an intermediate region.

2.4 Asymptotics of the probability density in the sticky limit

We pursue the method of matched asymptotic expansions to show that the leading-order probability density of $X^{\epsilon}_{t}$ (leading-order in the sense of a measure) is

[TABLE]

where $p(x,t)$ solves

[TABLE]

with a given initial condition $p(x,0)=\phi(x)$. We assume that $\phi$ is continuously differentiable on $[0,\infty)$ and that $\phi$ and $\phi^{\prime}$ are bounded on $[0,\infty)$; assumptions which guarantee that (6) is well-posed [79, Theorem 1]. Equations (6), (7) are the main results of this section, and describe the evolution of probability of a SBM [27].

To derive these equations we adopt a formal approach to highlight the main ideas, but expect that the argument could be turned into a more rigorous proof of weak convergence by the interested reader.

Let us make the ansatz that the solution to (4) away from the origin ($x\geq\epsilon$), or outer solution, and near the origin ($x\leq\epsilon$), or inner solution, have the asymptotic expansions, respectively,

[TABLE]

where $p_{0}(x,t)$ and $q_{0}(x,t)$ are the leading order terms in the expansions, and $p_{1}\ll p_{0}$, $q_{1}\ll q_{0}$ as $\epsilon\to 0$. We do not assume any particular scaling for these terms, nor the leading-order solutions; it will turn out that $p_{0}\sim O(1)$ but $q_{0}\sim O(\epsilon^{-1})$, because of Assumption 2.2 (A3).

By Assumption 2.2 (A1), the outer solution $p_{0}(x,t)$ satisfies the following linear PDE to leading order as $\epsilon\to 0$:

[TABLE]

with $p_{0}(x,0)=\phi(x)$. Near the origin, the density changes rapidly so it is convenient to change variables to $X=x/\epsilon$. Keeping the leading-order terms gives

[TABLE]

with the reflecting boundary condition $\partial_{X}U^{\epsilon}(0,t)q_{0}(0,t)+\partial_{X}q_{0}(0,t)=0$. We have used Assumption 2.2 (A2) and the timescale it implies to argue that $\partial_{t}q_{0}$ is lower order than the terms retained in (10).

Equation (10) is the equation for the stationary density of a particle diffusing in a potential $U^{\epsilon}$ with a reflecting boundary condition at the origin. The solution is

[TABLE]

where $a(t)$ is some unknown function of time, to be determined by matching to the outer solution.

Now we match the outer solution $p_{0}(x,t)$ and the inner solution $q_{0}(\epsilon X,t)$ at $x=\epsilon$, $X=1$. Unlike the traditional method of matched asymptotic expansions, we do not match in an overlap region, but rather at a single point, which is possible because the perturbation is not singular (the diffusion terms do not change with $\epsilon$.) We match using two conditions: one, probability is continuous, and two, probability is conserved. The first condition requires that

[TABLE]

to leading order in $\epsilon$, where we used that $e^{-U^{\epsilon}(\epsilon)}\sim 1$, by Assumption 2.2 (A1), and $p_{0}(\epsilon,t)\sim p_{0}(0,t)$. The condition that probability be conserved requires that

[TABLE]

to leading order in $\epsilon$. We moved the time derivative into the integrals, substituted for $\partial_{t}p_{0}$ using (9), and substituted $\kappa$ using Assumption 2.2 (A3).

Putting these results together gives a boundary condition for the outer solution $p_{0}$ at the origin, which can be written in two ways:

[TABLE]

Combining with (9) and removing the subscript on $p$ gives (7).

To obtain (6), notice that by continuity of probability, the leading-order density is $p_{0}(x,t)e^{-U^{\epsilon}(x)}$. Assumption 2.2 implies that $e^{-U^{\epsilon}(x)}$ converges weakly to $1+\kappa\delta(x)$ as $\epsilon\to 0$, giving a limiting density of (6).

2.5 Fokker-Planck equation

We make a few remarks concerning the limiting dynamics (7) and their relationship to a sticky Fokker-Planck equation.

Notice that, as expected, the limiting probability density (6) is singular with respect to Lebesgue measure on $\mathbb{R}^{+}$. Its structure shows that an SBM can spend finite time on any interval in its domain, as well as at the origin, $\{0\}$, sets with intrinsically different dimensions. This is an unusual property for diffusion processes; it does not occur for processes with the more typically-studied Dirichlet, Neumann, or Robin boundary conditions.

It turns out that the measure $1+\kappa\delta(x)$ in (6) is infinitesimally invariant, and is proportional to the invariant measure if the SBM is confined to a compact space, as the following lemma shows.

Lemma 2.5**.**

Suppose the function $p$ in (6), (7) additionally satisfies a reflecting boundary condition $p_{x}=0$ at $x=L>0$. Then the unique invariant probability measure $\pi$ is

[TABLE]

This is also the measure obtained as the weak limit of the stationary densities for (4).

Proof 2.6**.**

*Solve for the steady-state solution of (7) to obtain $p(x,t)=c$ for some constant $c\in\mathbb{R}$, and use (6) to obtain (13). *

It is sometimes helpful to work with a formulation of (7) directly in terms of the probability density. This is possible by first writing the density $\rho$ as a sum of densities on the different manifolds in its support, as

[TABLE]

where $\mu_{0}(dx)$ is the delta-function measure on $\{0\}$, and $\mu_{1}(dx)$ is the Lebesgue measure on $(0,\infty)$. We must impose a ‘continuity condition’

[TABLE]

to be consistent with the asymptotic derivation. Calculating $\partial_{t}\rho$ using (7) and (6) gives the system of equations

[TABLE]

System (16) and the continuity condition (15) can be interpreted as the Fokker-Planck equation for the evolution of the probability density $\rho$. If we write the system formally as $\partial_{t}\rho=\mathcal{L}^{*}\rho$ for some linear operator $\mathcal{L}^{*}$, then (16) shows that $\mathcal{L}^{*}$ should, formally at least, be interpreted as being a different partial differential operator, for the densities on each different manifold in the support of the process. Indeed, it turns out that for higher-dimensional sticky diffusions one can impose different dynamics in the interior of a domain and on the boundary, and these dynamics don’t have to bear any relation to each other [44].

The system (16) gives a more physical interpretation of the boundary condition in (7). It shows the boundary condition simply balances fluxes: the rate of change of the probability at the origin, $\partial_{t}\rho_{t}|_{x=0}=\partial_{t}(\kappa p)|_{x=0}$, equals the flux of probability which leaves the open interval $(0,\infty)$ on the left, $\partial_{x}\rho_{1}|_{x=0}=\partial_{x}p|_{x=0}$. Here again we see how the sticky boundary condition interpolates between a reflecting condition, when $\kappa=0$ and the condition is $\partial_{x}p|_{x=0}$, and an absorbing one, when $\kappa\to\infty$, and the condition approaches $\partial_{t}p|_{x=0}=0$.

2.6 Generator

We now consider how to obtain the generator $\mathcal{L}$ of a SBM, starting from the Fokker-Planck equation (16). The generator forms the basis for our numerical method, and is in fact the more fundamental quantity describing a Markov process. Recall that the generator is the operator $\mathcal{L}$ which is the formal adjoint of $\mathcal{L}^{*}$, i.e. $\langle f,\mathcal{L}^{*}\rho\rangle=\langle\mathcal{L}f,\rho\rangle$, for all densities $\rho$ of the form (6), and all test functions $f\in C^{2}_{c}([0,\infty))$ that satisfy an appropriate boundary condition at zero (to be determined.) Because properly defining $\mathcal{L}^{*}$ is somewhat subtle due to the singularities at the origin, we find it more transparent to work with a weakly equivalent formulation, which asks that $\mathcal{L}$ satisfy

[TABLE]

To proceed, compute the left-hand side of (17) as

[TABLE]

In the first step we substituted $\partial_{xx}p$ for $\partial_{t}p$ and integrated by parts, assuming a decay condition at $x=\infty$, and in the second step we used the boundary condition on $p$ at $x=0$ and rewrote the integral in terms of $\rho$.

We see that (17) is satisfied if we choose the generator and its associated boundary conditions to be

[TABLE]

For a SBM, the generator and its formal adjoint are equal.

Notice that if in (17) we had replaced $\rho$ with $p$, the function solving (7), we would not have found appropriate boundary conditions for $f$. This shows the importance of interpreting the weak formulation of (7) in the correct function space.

From the generator we obtain another way to verify that (13) is an infinitesimally invariant measure, by showing $\langle\mathcal{L}f,\pi\rangle=0$ for all test functions $f$ satisfying the appropriate boundary conditions. Calculations very similar to the above show this equation holds.

2.6.1 Example: MFPT of a sticky Brownian motion

As a simple application of (18), we directly calculate the MFPT of an SBM.

Lemma 2.7**.**

*The MFPT of an SBM out of $[0,\ell]$ with initial condition $x\in[0,\ell]$ is $\tau(x)=-\kappa x-(1/2)x^{2}+\kappa\ell+(1/2)\ell^{2}$. In particular, $\tau(0)=\kappa\ell+(1/2)\ell^{2}$. *

Proof 2.8**.**

The MFPT of an SBM satisfies the boundary value problem $\mathcal{L}\tau(x)=-1$, $\tau(l)=0$, plus any other boundary conditions associated with $\mathcal{L}$. Specifically,

[TABLE]

*The solution is $\tau(x)=-(1/2)x^{2}+c_{1}x+c_{2}$ where $c_{1}$ and $c_{2}$ are constants determined by imposing the boundary conditions. *

Note that $\tau(0)=\lim_{\epsilon\to 0}\tau^{\epsilon}(0)$, the limit of the MFPT for the reflecting diffusions, as in Lemma 2.3.

2.6.2 Example: transition rates between sticky points

As another application of (18), we consider a Brownian motion on a line segment with sticky endpoints, and calculate the transition rates between the endpoints. Such a setup is a good model for the transition paths between clusters of mesoscale particles, and quantitatively predicts their transition rates [78]. Suppose the endpoints have sticky parameters $\kappa_{1},\kappa_{2}$ and the line has length $L$, as shown below.

\kappa_{1}$$\kappa_{2}$$L[math]

The generator for this problem is

[TABLE]

The sign of $\partial_{x}f$ is reversed for the boundary condition at $L$, because the probability flux is in the opposite direction from the flux at [math].

A framework for calculating transition rates between disjoint sets $A$, $B$ exactly is given by Transition Path Theory [19]. These transition rates can be determined from empirical averages by using the following limit relations

[TABLE]

where $T$ is the time of observation, $N_{T}$ is the total number of transitions observed from $A$ to $B$ in time $T$, and $T_{A},T_{B}$ are the total times during which the process last hit $A,B$ respectively; they satisfy $T_{A}+T_{B}=T$. We apply Transition Path Theory to the sets $A=\{0\}$ and $B=\{L\}$ to compute the quantities above, and refer the reader to [19] for more justification of these calculations.

Lemma 2.9**.**

Given a process with generator $\mathcal{L}_{L}$ defined in (19), and let $A=\{0\}$, $B=\{L\}$. Then the stationary distribution is $\pi(x)=\frac{1}{L+\kappa_{1}+\kappa_{2}}(1+\kappa_{1}\delta(x)+\kappa_{2}\delta(x-L))$, and the transition rates are

[TABLE]

It is worth noting that the transition rates depend inversely on $L$, the total distance that the process must diffuse. This is a contrast to theories of transition rates derived from the Arrhenius formula or Transition State Theory, for which the transition rates depend only on local properties of critical points, such as the energy difference between critical points and their curvatures [32]. In the problem above, the energy differences are contained solely in $\kappa_{1},\kappa_{2}$.

Proof 2.10**.**

One can verify that $\langle\mathcal{L}_{L}f,\pi\rangle=0$ for all twice-differentiable functions $f$ satisfying the boundary conditions, so $\pi$ given above is the stationary distribution. Now we turn to calculating the transition rates. We first calculate the committor function $q(x)$, a function that gives the probability of hitting $B$ first before $A$ starting from $x$. It solves the boundary value problem

[TABLE]

*whose solution is $q(x)=x/L$. Next we calculate the overall rate of transition $\nu=\lim_{T\to\infty}\frac{N_{T}}{T}$, as $\nu=\pi(x)\partial_{x}q(x)$ for any $x\in(0,L)$, giving $\nu=L^{-1}(L+\kappa_{1}+\kappa_{2})^{-1}$. Then we calculate the “reactive probabilities” $\rho_{A}=\lim_{T\to\infty}\frac{T_{A}}{T}$, $\rho_{B}=\lim_{T\to\infty}\frac{T_{B}}{T}$. These are computed as $\rho_{A}=\int_{0}^{L}\pi(x)(1-q(x))dx$, $\rho_{B}=\int_{0}^{L}\pi(x)q(x)dx$, giving $\rho_{A}=(\kappa_{1}+L/2)/(L+\kappa_{1}+\kappa_{2})$, $\rho_{B}=(\kappa_{2}+L/2)/(L+\kappa_{1}+\kappa_{2})$. Finally, the reaction rates are computed as $k_{A\to B}=\nu/\rho_{A}$, $k_{B\to A}=\nu/\rho_{B}$, giving the result above. *

3 Numerical method to simulate a sticky Brownian motion

It is not obvious from either the Fokker-Planck equation or the generator how one should simulate trajectories of a sticky Brownian motion. One option is to return to the derivation as a limit as $\epsilon\to 0$ of $X^{\epsilon}_{t}$, and choose a small $\epsilon$ and simulate (1). However, our earlier calculations (Figure 2) showed that for small $\epsilon$ one must use a very small timestep; for large $\epsilon$ the solution will be inaccurate.

Therefore, we turn to an entirely different method, based on constructing a continuous-time Markov chain whose generator $Q$ approximates the generator $\mathcal{L}$ of an SBM. Specifically, we spatially discretize the infinitesimal generator and boundary conditions of the SBM, for example using a finite difference approximation, to obtain the generator of a Markov jump process on the set of discretization points. This Markov jump process may be simulated using a simple Monte Carlo method known variously as the Stochastic Simulation Algorithm, kinetic Monte Carlo, the Doob-Gillespie algorithm or the Gillespie algorithm [16, 17, 33, 34].

These approximations of sticky diffusions go back to the Markov Chain Approximation Method (MCAM) invented by Harold Kushner in the 1970s to approximate optimally controlled diffusion processes [51, 61, 62, 52, 53, 54, 55, 60, 56, 57, 58, 59]. However, because of their interest in stochastic control problems, these works mainly focus on numerical solutions with gridded state spaces and use numerical linear algebra to construct an approximation. In the physics literature, a Monte-Carlo method was developed to construct the approximation [21, 89, 73, 74, 63], an approach which seems to go back to at least [21]. More recently, a new “gridless” framework was introduced for constructing Markov jump process approximations for diffusions [9], which allows the domain of the diffusion process to be unbounded, does not require that the diffusion process is symmetric, and does not assume that the infinitesimal covariance matrix of the diffusion process is diagonally dominant. This generalization allows the jump size of the numerical solution to be uniformly bounded, which makes it easier to numerically treat boundary conditions. Very recently, these ideas have been extended to SPDE problems [8].

3.1 Numerical Algorithm

To apply this approach, let us discretize the interval $[0,\infty)$ into a set of grid points $\mathcal{R}^{h}=\{0,h,2h,\ldots\}$, where $h$ is the spacing between neighboring grid points. Let $f\in C^{4}_{b}(\mathbb{R})$ and let $f_{0},f_{1},\ldots$ be the values of the function at the grid points, i.e. $f_{k}=f(kh)$ as illustrated below.

$-h$[math]h$$2h$$f_{-1}$$f_{0}$$f_{1}$$f_{2}$$\dotsc

The white dot indicates a ghost grid point, which we recruit in our construction.

At each interior grid point, the generator $\mathcal{L}f$ (see (18)) can be approximated by

[TABLE]

This approximation is a second-order centered finite difference approximation to the second derivative operator $\partial_{xx}$. At the boundary grid point, $k=0$, we do not know the value of the ‘ghost’ point $f_{-1}$ that is needed in (20), so we solve for it using the boundary condition in (18). The discretized boundary condition is

[TABLE]

where again we used a second-order centered finite difference scheme to evaluate $\partial_{x}f$. Using this equation to eliminate $f_{-1}$ in (20) with $k=0$ gives an approximation to the generator at the boundary as

[TABLE]

This local error estimate suggests that this approximation is globally second-order accurate in $h$ when $\kappa>0$. When $\kappa=0$, which corresponds to a reflecting Brownian motion, the remainder term in (21) may lead one to surmise that the approximation is only first-order spatially accurate. However, it turns out that even in this case, the global error is still second-order in $h$, because roughly speaking, the local error at the boundary is also proportional to the mean occupation time of a reflecting Brownian motion near the boundary, which is $O(h)$.

Now, we construct a continuous-time Markov chain $Y_{t}$ with state space $\mathcal{R}^{h}$, whose generator $Q$ is the discrete approximation to $\mathcal{L}$. That is, we set

[TABLE]

The infinite matrix associated to the generator $Q$ has nonzero entries $Q_{kj}=1/h^{2}$ if $j=k\pm 1$ and $k\geq 1$, $Q_{kk}=-2/h^{2}$ if $k\geq 1$, $Q_{01}=2/(h^{2}+2\kappa h)$, $Q_{00}=-2/(h^{2}+2\kappa h)$. One can verify that $Q$ is the generator of a continuous-time Markov chain, since the coefficients $Q_{kj}$ with $k\neq j$ are nonnegative, and $\sum_{j}Q_{kj}=0$. Note that when $\kappa=0$ in (22) one obtains the generator of a reflecting random walk.

Realizations of this sticky random walk $Y_{t}$ can be simulated exactly using a simple Monte Carlo method [16, 17, 33, 34]. Suppose $Y_{t_{0}}=kh$ at time $t_{0}$. The process is updated in two steps.

• •

Pick a state $\nu\in\{(k+1)h,(k-1)h\}$ with probabilities proportional to $Q_{k,k+1}$, $Q_{k,k-1}$ respectively. Specifically, for $k=1,2,\ldots$, we jump to the left or right states with equal probability, and for $k=0$, we always jump to $k=1$.

• •

Choose a random time $\tau\sim\mbox{Exp}(-Q_{kk})$, where $\mbox{Exp}(\lambda)$ is the distribution of an exponential random variable with mean $1/\lambda$. Set $Y_{t_{0}+t}=kh$ for $0<t<\tau$ and set $Y_{t_{0}+\tau}=\nu$. That is, the process jumps to a new state $\nu$ after a random time $\tau$.

To sum up, the process $Y_{t}$ is a continuous-time random walk, a “sticky random walk” (SRW), with mean waiting time $h^{2}/2$ at interior grid points and $\kappa h+h^{2}/2$ at the sticky boundary. See Appendix A for Matlab code implementing this method.

Realizations of $Y_{t}$ for $\kappa=1$ and different values of $h$ are plotted in Figure 4. Away from 0, the process is a random walk with random waiting times, which looks increasingly like a Brownian motion as $h\to 0$. When the process hits 0, however, it spends much longer there on average than it does at interior points. Overall the process looks like a Brownian motion which has been “slowed down” near 0, an observation that can be made rigorous as we discuss in Section 4.

3.2 Properties of the numerical solution

This difference in the scaling with $h$ of the waiting times at boundary and interior points, nicely illustrates the unusual behavior of an SBM at the origin and gives insight into why it is sticky there. The mean waiting time at interior points is $h^{2}/2$, the usual scaling for diffusion motion, which requires that $\mathbb{E}(\Delta Y)^{2}/\tau\to cst$, where $\Delta Y=\pm h$ is the jump at each step of the algorithm. The mean waiting time at the origin is $\kappa h+h^{2}/2=\kappa h+o(h)$ – longer by a factor of about $h$, so the scaling is ballistic at the origin. The amount of time the process spends at the origin is the same as it would spend at about $\kappa/h$ interior grid points, which is the number one would use to discretize an interval of length $\kappa$. Therefore, we expect we expect this SWR to spend a finite, non-zero amount of time at the origin as $h\to 0$, an amount of time which increases in proportion to $\kappa$.

Interestingly, although the SRW spends finite time at the origin as $h\to 0$, in the limit it never spends a whole interval of time at the origin. This is evident from the numerical solution, since the waiting time at the origin approaches zero as $h\to 0$. This implies that every time an SBM hits the origin, it leaves right away, just as does a Brownian motion; yet somehow the total measure of the points at which it equals zero is positive.

We have seen the ballistic scaling of times near the origin before, in Lemmas 2.3 and 2.7. Indeed, the mean holding time of the sticky random walk $Y_{t}$ at [math] is identically the MFPT of a SBM with $\ell=h$ (see Lemma 2.7.) The additive property of the MFPT then implies the following lemma.

Lemma 3.1**.**

*Given $h>0$ and $\ell\in\mathcal{R}_{h}$, the MFPT of $Y_{t}$ out of $[0,\ell]$ with initial condition $x\in[0,\ell]\cap\mathcal{R}_{h}$ is $\tau(x)=-\kappa x-(1/2)x^{2}+\kappa\ell+(1/2)\ell^{2}$. In particular, $\tau(0)=\kappa\ell+(1/2)\ell^{2}$. *

Another useful property of $Y_{t}$ is that once a simulation has been performed for some fixed $\kappa$, one can obtain a trajectory for any other value $\kappa^{\prime}\neq\kappa$ simply by scaling the holding time of the jumps which leave the origin, since the probabilities of jumping to each state do not change; see Figure 5 for an illustration. This is a useful property as it allows one to investigate multiple values of $\kappa$ with one single simulation.

3.3 Accuracy and efficiency of the numerical method

We compare the accuracy and efficiency of this numerical method for simulating a SBM, to the accuracy and efficiency of a method which simulates (1) directly using a symmetrized Euler-Maruyama (EM) approximation (5) and a Morse potential for $U^{\epsilon}(x)$. Figure 6 shows that 1% accuracy for the sticky random walk requires a mean time step of about 0.1. Contrast this to the EM scheme, which, for an accuracy of about 2%, requires a potential width no more than about $118^{-1}$ and a time step of about $1.25\times 10^{-6}$. In this example, the time step of the sticky random walk method is about five orders of magnitude larger, to obtain a comparable level of accuracy.

In applications, one may with to use a SBM to model a potential with a small, but not infinitesimal range. How much more efficient is the sticky random walk in capturing the statistics of (1) for finite $\epsilon$, and how accurately does it do so? To test this we ran a high-resolution simulation of (1) ($\delta t=3.91\times 10^{-6}$) using a Morse potential (Example 2.1) with sticky parameter $\kappa=30$, and range parameters $a\approx 30$ (well depth $D_{e}=7.22$), characteristic of certain depletion interactions [72, 78], and $a\approx 100$ ($D_{e}=8.5$), characteristic of certain DNA-mediated interactions [90]. We computed the statistics $\mathbb{E}X_{1}^{\epsilon}$, $\mathbb{P}(X_{1}^{\epsilon}<0.15)$, and compared them to statistics estimated using the sticky random walk method with $\kappa=30$. Figure 7 shows that the sticky random walk method gives virtually the same estimates for all time steps; in particular the largest mean time step used, $\delta t\approx 10^{-2}$, is as good as any smaller timestep. In contrast, statistics computed from a direct simulation of (1) are nowhere near their converged values until a timestep of $\delta t=3\times 10^{-5},2\times 10^{-6}$ for $a=30,100$ respectively, more than 2 orders of magnitude smaller than for the sticky random walk.

The sticky random walk method will always contain a non-zero relative error no matter what the timestep, since it is only an approximation to the true dynamics. How large is this error? For the mean, $\mathbb{E}X_{1}^{\epsilon}$, the relative error between the sticky random walk and the high-resolution simulation is $40\%,13\%$ for $a=30,100$ respectively, and for $\mathbb{P}(X_{1}^{\epsilon}<0.15)$, the relative errors are $6\%,2\%$. While the relative error in computing the mean may seem high, notice that the mean is sensitive to the range and shape of the potential, since we are computing this statistic at short times so most of the probability is concentrated in the boundary layer associated with the potential. In practice, experimental measurements of a quantity that depends so sensitively on the shape of a short-ranged potential are infeasible, because measurement noise is typically larger than the width of the potential [78]. In addition, the quantity of interest would usually be the particle diameter plus the mean distance, which would have a significantly smaller relative error, less than 1% for a 1$\mu$m particle. A coarser statistic such as $\mathbb{P}(X_{1}^{\epsilon}<0.15)$ is more practical to measure, and hence is a better test for the accuracy of the sticky random walk model; here the relative errors are much smaller.

3.4 PDEs that can be solved by running a sticky random walk

It should be clear from the calculations in section 3 that the generator $Q$ of $Y_{t}$ is a locally consistent approximation of the generator $\mathcal{L}$ of a SBM. How does this local consistency connect to more global results, for example about the accuracy with which $Y_{t}$ can reproduce statistics of a sticky Brownian motion? In this section we present new Feynman-Kac formulae to show that from trajectories of $Y_{t}$, one can numerically solve certain PDEs with Feller’s more general boundary condition to second-order accuracy in the space step.

Consider the heat equation with Feller’s boundary condition

[TABLE]

where $\phi:[0,\infty)\to\mathbb{R}$ and $p_{1},p_{2},p_{3}$ are parameters satisfying $p_{1}+p_{2}+p_{3}=1$ and $p_{1},p_{2},p_{3}\geq 0$ [27, 29, 45]. Special cases of this boundary condition include:

1) $p_{1}=1$:

Dirichlet b.c. corresponding to a stopped Brownian motion;

2) $p_{2}=1$:

Neumann b.c. corresponding to a reflecting Brownian motion;

3) $p_{3}=1$:

Wentzell b.c. corresponding to an absorbed Brownian motion;444An absorbed Brownian motion is one where the Brownian motion has a fixed point at the boundary, which is consistent with what happens as the sticky parameter in (7) becomes infinite; whereas, a stopped Brownian motion corresponds to a Brownian motion that is terminated at the boundary. and,

4) $p_{3}=0$:

Robin b.c. corresponding to an elastic Brownian motion.

The general case corresponds to a Brownian motion that exhibits stopping, stickiness, and reflection at various rates [45].

The following result shows how to numerically solve (23) when $p_{2}>0$ using the sticky random walk developed in §3. The case $p_{2}=0$ can be solved in a similar manner using a continuous-time random walk with absorbing boundary conditions at [math], but for notational brevity, we omit this case.

Theorem 3.2**.**

Assume $p_{2}>0$, set

[TABLE]

Suppose further that (23) has a solution $u(x,t)$ whose first four derivatives are continuous and bounded. Let $Y_{t}$ be a sticky random walk, i.e., a Markov jump process on the grid $\mathcal{R}^{h}$ with generator $Q$ defined by (22). At every grid point $x_{i}\in\mathcal{R}^{h}$, define the function

[TABLE]

where $\mathbb{E}_{x_{i}}$ denotes expectation conditional on $Y_{0}=x_{i}$. Then for all $T>0$ and for all $h>0$ sufficiently small, there exists $C(T)>0$ such that

[TABLE]

Although there are well-posedness and semigroup generation results for the heat equation with generalized, two-sided Feller boundary conditions in a bounded domain [24, 25], we are not aware of results for the heat equation on the half line with Feller boundary conditions that gives the required regularity results for $u(x,t)$ in Theorem 3.2.

Proof 3.3**.**

Our proof will proceed in two steps. First, we show that a second-order spatial discretization of the PDE (23) using finite differences is given by the solution $u^{h}(t)\in\mathcal{R}^{h}$ to the linear, infinite-dimensional system of ODEs

[TABLE]

with initial condition $u^{h}(0)=\left.\phi\right|_{\mathcal{R}^{h}}$. Together with the regularity of $u(x,t)$, this will imply that $u^{h}_{i}(t)$ and $u(x_{i},t)$ are close in an $\ell_{\infty}$ norm over $\mathcal{R}^{h}$. Second, we show that the solution $u^{h}_{i}(t)$ to (25) admits the stochastic representation (24).

Note that another way to write (25) is to define $u^{h}(x_{i},t)=u_{i}^{h}(t)$, and write

[TABLE]

where $Q$ is the generator of $Y_{t}$ defined in (22). The difference between (26) and the backward equation associated with $Y_{t}$ is the presence of the term $cu_{0}^{h}$, which causes $u(0,t)$ to decrease in time (recall $c<0$).

For the first step, i.e. obtaining the discretization of the PDE (28), notice that the ODEs (25) for $k\geq 1$ are obtained by replacing $\partial_{xx}$ in (23) by the centered finite difference given in (20); and the ODE at $k=0$ is obtained by incorporating a finite difference discretization of the boundary condition into a centered finite difference approximation of $\partial_{xx}$ at [math]. Specifically, given a function $f$ that satisfies the boundary condition in (23), use a ghost grid value $f_{-1}$ combined with a centered finite difference approximation to $f^{\prime}(0)$ and $f^{\prime\prime}(0)$ to obtain

[TABLE]

Using the above equation to eliminate $f_{-1}$ from (20) when $k=0$, gives an approximation to $\partial_{xx}$ at the boundary as

[TABLE]

Substituting $p_{3}=\kappa p_{2}$ into the first term then gives the ODE in (25) for $k=0$.

For the second step, define

[TABLE]

Applying Itô’s formula for jump processes (see e.g. [82]) we obtain that

[TABLE]

By (26), the right-hand side is zero. Evaluating this equation at $s=t$ and using the initial condition in (23) gives

[TABLE]

*as required. *

One can similarly solve Dirichlet-Poisson problems. For example, consider the boundary value problem

[TABLE]

where $\phi:[0,\ell]\to\mathbb{R}$ is a given function and $p_{1},p_{2},p_{3}$ are parameters satisfying $p_{1}+p_{2}+p_{3}=1$ and $p_{1},p_{2},p_{3}\geq 0$. Then we have the following result.

Theorem 3.4**.**

Assume $p_{2}>0$, set $\kappa=p_{3}/p_{2}$ and $c=-2p_{1}/(hp_{2}+2p_{3})$; and suppose that (28) has a solution $u(x)$ whose first four derivatives are continuous and bounded. Let $Y_{t}$ be a sticky random walk, i.e., a Markov jump process on the grid $\mathcal{R}^{h}$ with generator $Q$ defined by (22). For every grid point $x_{i}\in\mathcal{R}^{h}\cap[0,\ell]$, define the function

[TABLE]

where $\mathbb{E}_{x_{i}}$ denotes expectation conditional on $Y_{0}=x_{i}$ and $\tau^{h}$ is the random stopping time $\tau^{h}=\inf\{t\geq 0\mid Y_{t}=\ell\}$. Then for all $h>0$ sufficiently small, there exists $C>0$ such that

[TABLE]

A proof of this theorem is like the proof of Theorem 3.2 and therefore omitted.

4 Other characterizations of sticky Brownian motion

Interestingly, sticky boundary conditions were discovered in the pure mathematics literature, well before they were used in applied math or physics. This section summarizes some of the other ways of characterizing SBM that have been developed in the mathematical literature, and connects them to our numerical and asymptotic approach. A nice history focusing especially on Feller’s contribution to the development is given in [80]. To more easily make a connection to the existing literature, we consider in this section a traditional sticky Brownian motion, whose generator is $\mathcal{L}^{0}f=\frac{1}{2}\partial_{xx}f$ with boundary condition $\partial_{x}f|_{x=0}=\kappa\partial_{xx}f|_{x=0}$. The numerical solution of this process is the same as that given in section 3, but with numerical generator $Q^{0}=Q/2$, where $Q$ is the generator for a root-2 sticky Brownian motion, defined in (22).

Sticky diffusion processes were discovered by Feller in the 1950s, who sought to identify the most general boundary conditions for the generator of a one-dimensional diffusion [27, 29]. For a process on $[0,\infty)$ with continuous sample paths that behaves like a Brownian motion in $(0,\infty)$, the possible boundary conditions have the form [27, 45, 80]

[TABLE]

with $p_{1},p_{2},p_{3}\geq 0$, and $p_{1}+p_{2}+p_{3}=1$ (c.f. (23)). Feller said of his result [27]:

This is the second instance of a concrete problem with physical significance where physical intuition failed, but our abstract methods provided a clue.

Feller was referring to problems in genetics, which have certain singularities in the diffusion coefficients when a population size hits zero, but his remark could equally well apply to the sticky boundary condition $f^{\prime}(0)=\frac{p_{3}}{p_{2}}f^{\prime\prime}(0)$; physical intuition does not immediately show why this should be associated with stickiness, and a singularity in the transition probabilities.

Subsequently Itô & McKean showed how to construct an SBM from a time change of a reflecting Brownian motion [45]. Let $B_{t}^{+}$ be a standard reflecting Brownian motion, which can be constructed from a standard Brownian motion $B_{t}$ as $B^{+}_{t}=|B_{t}|$. Let

[TABLE]

be the local time accumulated at [math] by $B^{+}$ over $[0,t]$. Define

[TABLE]

Since $l^{0}_{t}(B^{+})$ is continuous and nondecreasing, $A_{t}$ is continuous and strictly increasing, so it has an inverse $T(t)=A^{-1}(t)$. Define

[TABLE]

Then $B^{*}_{t}$ is an SBM [45].

The process $B^{*}_{t}$ is simply a time change of a reflecting Brownian motion, making it run on clock $T(t)$ instead of clock $t$. When $B^{*}_{t}\neq 0$, then $dA/dt=dT/dt=1$ so the process runs at the same rate as $B^{+}_{t}$. When $B^{*}_{t}=0$, then $A(t)$ increases faster than $t$ because of the local time, so the clock $T(t)$ runs more slowly and $B^{*}_{t}$ slows down at [math].

To see why it slows down in a way that changes its measure at zero, let’s compute the MFPT for $B^{*}_{t}$ to leave the interval $[0,h]$, starting at 0. Let $\tau=\inf\{t>0:B_{t}^{+}>h\}$ be the first exit time of $B_{t}^{+}$ from $[0,h]$, so that $\mathbb{E}(\tau)=h^{2}$ and the MFPT of $B_{t}^{*}$ is given by

[TABLE]

By a Feynman-Kac formula, the expected value $u(x)=\mathbb{E}_{x}(\int_{0}^{\tau}1(B_{s}^{+}<\epsilon)ds)$ can be calculated by solving the boundary value problem

[TABLE]

whose local solution at zero is $u(0)=2\epsilon h-\epsilon^{2}$. Hence, $\mathbb{E}(\ell_{\tau}^{0}(B))=h$ and $\mathbb{E}(A_{\tau})=h^{2}+2\kappa h$, which is twice the mean waiting time we derived for a root-2 SBM in Lemma 2.7, and in particular, this waiting time is not diffusive like the mean waiting time of $B_{t}^{+}$.

A different discrete approximation to a SBM was derived by Amir [1], who showed how to obtain the process as a limit of random walks. Recall that one can construct a standard Brownian motion starting from a random walk by using Donsker’s Theorem [15]. Let

[TABLE]

be a rescaled random walk, where $\{X_{i}\}_{i=1}^{\infty}$ is a sequence of independent, identically distributed random variables with $P(X_{i}=1)=P(X_{i}=-1)=\frac{1}{2}$. As $\Delta x\to 0$ with $\Delta t=(\Delta x)^{2}$, the process $S_{\Delta x}(t)$ converges in distribution to a Brownian motion $B_{t}$ [15, 26].

Amir [1] showed that one can construct a sticky (nonreflecting) Brownian motion, by modifying $S_{\Delta x}$ as follows: every time $S_{\Delta x}(t)$ hits zero, wait there for a time interval of length $\sqrt{\Delta t}$, instead of $\Delta t$. For example, one can let $\Delta x=2^{-n}$, $\Delta t=2^{-2n}$ for some integer $n$, and then the process must wait $2^{n}$ timesteps each time it hits zero. As $n\to\infty$, the modified process $S^{*}_{n}(t)$ converges to a sticky nonreflecting Brownian motion. A sticky reflecting Brownian motion is then the limit of $|S^{*}_{n}(t)|$.

Finally, we remark that an SBM $B^{*}_{t}$ solves the following SDE, which should be interpreted in an integrated sense (see [44, 22], and references therein):

[TABLE]

This equation has a unique weak solution, but no strong solution [22]. That it corresponds to the sticky boundary condition can be seen by noting that555This perspective is due to R. Varadhan, personal communication. the generator away from the origin is $\mathcal{A}f=\frac{1}{2}\partial_{xx}f$, and at the origin it is $\mathcal{A}f=\frac{1}{2\kappa}\partial_{x}f$, so to be consistent the function $f$ must satisfy the sticky boundary condition $\frac{1}{2}\partial_{xx}f=\frac{1}{2\kappa}\partial_{x}f$.

It may seem quite remarkable that changing the diffusion coefficient at a single point, the origin, can have such a dramatic effect on the process. To see why this is so, consider widening the region near the origin by some amount $\epsilon\ll 1$, approximating the drift as constant in this interval, and estimating the increment $\Delta B_{t}^{*}$ over a small time interval $\Delta t$. If $B^{*}_{t}<\epsilon$, the increment is $\Delta B^{*}_{t}\approx\frac{1}{2\kappa}\Delta t$, so it takes the process a time of about $\Delta t\approx 2\kappa\epsilon$ to leave the interval $[0,\epsilon)$. If $B^{*}_{t}>\epsilon$, the increment has magnitude $|\Delta B^{*}_{t}|\approx|\Delta B_{t}|\approx\sqrt{\Delta t}$, so it takes the process a time of about $\Delta t\approx\epsilon^{2}$ to travel the length of any other interval of length $\epsilon$. We recover the same scalings as for the time-changed formula (31), showing that the slowdown near the origin has a singular effect on its transition densities. Even though the diffusion coefficient changes at only one point, the difference between $dt$ and $dB_{t}$ at that point gives rise to the singular change in timescales.

5 Conclusion and Outlook

We considered a reflecting Brownian motion on a half-line with a deep but short-ranged potential energy near the origin. As the potential becomes narrower and deeper, the process approaches a “sticky” Brownian motion, which has finite probability to be found in any interval on the half-line, as well as finite probability to be found exactly at $\{0\}$; counterintuitively it never spends an interval of time at 0. The process is characterized by a non-classical boundary condition for its generator, involving second derivatives. We showed that simulating trajectories of a process that is close to sticky using a traditional Euler-Maruyama discretization of its SDE requires very small timesteps. This motivated us to introduce an alternative method, based on discretizing space first and constructing a continuous-time Markov chain on the set of discretization points, which allows time steps at least 2 orders of magnitude larger in parameter regimes of physical relevance (at the expense of making a small error in estimating the transition probability.) The method results in a random walk on the nonnegative integers with random holding times, with a larger holding time at [math]. The holding times give insight into why the process has a singular probability density at zero, since the holding time scales ballistically with step size at the origin, rather than diffusively.

Our motivation came from studying systems of mesoscale particles (diameters $\approx 100$nm-$10\mu$m), which have attractive interactions that are very short-ranged compared to their diameters. For such particles, the structures and dynamics are often more effectively studied by considering the system in the sticky limit [42, 78, 41]. Our goal is to eventually create numerical methods to simulate such particle systems directly in the sticky limit, for which we hope to see a similarly large gain in efficiency.

Achieving this goal will require building on the ideas in this paper to handle higher-dimensional sticky diffusions with more complicated boundary conditions. There are several steps involved. The first is to create methods to handle $d$-dimensional diffusions that are sticky on a half-space of dimension $d-1$. The new ingredient here is that the diffusion can move directly along the boundary, with dynamics that are different from those in the interior. This will be feasible by discretizing the generators in the domain and on its boundary, as we have done in this paper. We expect such a method to be efficient even when $d$ is too large to solve PDEs numerically, since the generator must only be discretized locally, using two points per dimension, and not globally over the whole domain [9].

A second step is to adapt these methods to work on manifolds, since a collection of interacting sticky particles performs a diffusion on the manifold corresponding to certain fixed distance constraints [42]. It may be a challenge to retain the method’s second-order accuracy when on a manifold, however it may not be necessary to do so, since the sticky limit is an approximation anyways. Other methods have considered how to sample probability densities directly on manifolds and have shown themselves to be significantly more efficient than using short-range forces to keep a process near a manifold (e.g. [13, 64, 92].)

A third step, and perhaps the most challenging, will be to adapt the method to processes which are sticky on even lower-dimensional “corners”. For example, a $d$-dimensional process may stick to a $d-1$-dimensional boundary, from which it may stick to a $d-2$-dimensional boundary, and so on. Physically this could correspond to a system of particles forming one bond, then two, then three, and so on. The number of intersecting boundaries increases as one moves down in dimension, and dealing with the whole collection of boundaries together may require further approximation.

We have focused in this paper exclusively on one-dimensional sticky Brownian motion, in order to build intuition into this unusual process, and to introduce it to the applied math community. Along the way we discussed some of the connections to the probabilistic approaches to studying sticky diffusions. We hope these ideas will be useful for other researchers studying sticky diffusions in the myriad of contexts in which they may arise, from biology to materials science to finance to operations research, and other applications yet to be imagined.

Appendix A Listings

In the MATLAB file SEM.m in Listing LABEL:list:SEM, we apply the symmetrized Euler-Maruyama scheme to the SDE (1) with $U^{\epsilon}(x)$ being a Morse potential energy function (see Example 2.1) with parameters $D_{e}=8.5$, $a=\sqrt{\pi}e^{D_{e}}/(\kappa\sqrt{D_{e}})$, and $x_{0}=1/a$; and the sticky parameter set at $\kappa=30$. A single realization of a discretized Brownian motion is produced over $[0,1]$ with $\delta t=10^{-6}$.

This realization of Brownian motion is used to drive symmetrized Euler-Maruyama operated at two time step sizes: $\delta t$ and $R\delta t$ where $R=16$ [39]. We set the seed for MATLAB’s random number generator, arbitrarily, to be 88 using the command rng. The plotted paths illustrate that the time step size of symmetrized Euler-Maruyama has to be sufficiently small in order for the method to accurately represent the strong, short-ranged Morse potential force.

Listing LABEL:list:SRW displays the MATLAB file SRW.m which produces a realization of the sticky random walk over the time interval $[0,T]$ with initial condition $x_{0}=0.2$, sticky parameter $\kappa=30$, and spatial step size $h=0.01$. The simulation is terminated after the time update first exceeds $T$. We set the seed for MATLAB’s random number generator, arbitrarily, to be 999 with the command rng. A single sample paths is plotted as a stairstep graph using the command stairs.

Acknowledgments

We wish to acknowledge Robert Kohn, Eric Vanden-Eijnden, and Srinivasa Varadhan for useful discussions that helped to improve this paper.