An encounter-based approach for restricted diffusion with a gradient drift

We first recall two conventional descriptions of restricted diffusion in a bounded domain $\Omega\subset{\mathbb{R}}^{d}$ with a smooth boundary $\partial\Omega$. On one hand, it can be described by the Skorokhod stochastic differential equation for the random position $\bm{X}_{t}=(X_{t}^{1},\ldots,X_{t}^{d})$ of a diffusing particle at time $t$:

where $\bm{x}_{0}$ is the starting point [9, 10, 11]. Here the infinitesimal displacement $d\bm{X}_{t}$ of the particle has three contributions: (i) the deterministic part with a drift $\bm{\mu}(\bm{x})$, (ii) random (thermal) fluctuations with the standard Gaussian noises $d\bm{W}_{t}$ whose amplitudes are determined by the diffusion coefficient $D$, and (iii) reflections from the boundary $\partial\Omega$ along the unit normal vector $-\bm{n}(\bm{x})$ (oriented inward the domain), whenever the particle hits the boundary. Intuitively, the last term can be understood as an infinitely strong but infinitely short-ranged force that repels the particle from the boundary back to the bulk (see a physical rational behind this force in SM.I of [8]). Such instant reflections are governed by the boundary local time $\ell_{t}$ – a non-decreasing stochastic process (with $\ell_{0}=0$) that increases at each encounter with the boundary (see Fig. 1). Curiously, the single Skorokhod equation determines simultaneously two tightly related stochastic processes: the position $\bm{X}_{t}$ and the boundary local time $\ell_{t}$. Once $\bm{X}_{t}$ is constructed, the boundary local time $\ell_{t}$ can be obtained by rescaling the residence time of the particle in the thin boundary layer of width $a$ (see [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25] and references therein):

where $\Theta(z)$ is the Heaviside step function ($\Theta(z)=1$ for $z>0$ and $0$ otherwise), and $|\bm{x}-\partial\Omega|$ denotes the distance between a bulk point $\bm{x}$ and the boundary $\partial\Omega$. From Eq. (2), one can see that the boundary local time $\ell_{t}$ is indeed a non-decreasing process that remains constant when $\bm{X}_{t}$ is the bulk, and increases only when $\bm{X}_{t}$ hits the boundary. Even though $\ell_{t}$ has units of length, we keep using the canonical term “boundary local time”. Note that $\ell_{t}/D$ has units of time per length so that its multiplication by the width $a$ of a thin boundary layer approximates the fraction of time that the particle spent in that layer up to time $t$.

Throughout this work, we focus on the common physical setting when the drift represents an external conservative force $\bm{F}(\bm{x})$, which can be written as the gradient of a potential $V(\bm{x})$:

where $\Phi(\bm{x})=V(\bm{x})/(k_{B}T)$ is the dimensionless potential and $\gamma=k_{B}T/D$ is the friction coefficient, with $k_{B}$ being the Boltzmann’s constant and $T$ the absolute temperature. We emphasize that the drift is time-independent while thermal fluctuations are isotropic and independent of both time and coordinates.

As said earlier, most physical works on restricted diffusion skip the stochastic differential equation (1) and deal directly with the probability density of the position $\bm{X}_{t}$ (also known as the propagator), $G_{q}(\bm{x},t|\bm{x}_{0})d\bm{x}={\mathbb{P}}_{\bm{x}_{0}}\{\bm{X}_{t}\in(\bm{x},\bm{x}+d\bm{x})\}$, which obeys the Fokker-Planck equation [2]

where

is the Fokker-Planck operator. Setting the probability flux density

the Fokker-Planck equation can be understood as the continuity equation expressing the probability conservation: $\partial_{t}G_{q}(\bm{x},t|\bm{x}_{0})=-(\nabla\cdot\bm{J}_{q})$.

The equation (4) has to be completed by the initial condition $G_{q}(\bm{x},0|\bm{x}_{0})=\delta(\bm{x}-\bm{x}_{0})$ stating that the particle has started from $\bm{x}_{0}$ at time $t=0$, and by an appropriate boundary condition that accounts for interactions with the boundary. When the boundary is inert and impermeable for the particle, the probability flux density in the normal direction to the boundary,

is zero at any boundary point $\bm{s}\in\partial\Omega$, with $\partial_{n}=(\bm{n}(\bm{x})\cdot\nabla)$ being the normal derivative oriented outwards the domain $\Omega$. In turn, a partially reactive boundary is often described by imposing that the probability flux density in the normal direction is proportional to $G_{q}(\bm{x},t|\bm{x}_{0})$,

with $\kappa$ (in units m/s) being the reactivity. This Robin-type boundary condition was introduced by Collins and Kimball in the context of chemical physics and was later investigated and employed by many authors [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41]. When $\kappa=0$, one retrieves the fully reflecting boundary, whereas in the opposite limit $\kappa\to\infty$, the left-hand side of this relation vanishes, and one gets the Dirichlet boundary condition $G(\bm{x},t|\bm{x}_{0})|_{\bm{x}\in\partial\Omega}=0$ for a perfectly reactive boundary. As discussed below, the Robin boundary condition (8) describes one surface reaction mechanism among many others.

It is well known that the propagator also satisfies the backward Fokker-Planck equation [2], which reads in our case as

where

is the adjoint Fokker-Planck operator. Since both $\bm{\mu}(\bm{x})$ and $D$ are time-independent, solutions of both Fokker-Planck equations depend on the difference between the terminal and starting times $t$ and $t_{0}=0$, that allowed us to replace the time derivative $-\partial_{t_{0}}$ by $\partial_{t}$ in Eq. (9).

The very definition of the propagator $G_{q}(\bm{x},t|\bm{x}_{0})$ indicates the deliberate choice of focusing on the position $\bm{X}_{t}$ and ignoring the boundary local time $\ell_{t}$ in the second description. As shown in [8] and discussed below, the inclusion of boundary encounters into the theoretical framework brings numerous advantages.

We start with restricted diffusion without any surface reaction and introduce the full propagator $P(\bm{x},\ell,t|\bm{x}_{0})$ as the joint probability density of $\bm{X}_{t}$ and $\ell_{t}$ in a bounded domain $\Omega$ with an impermeable boundary $\partial\Omega$:

In this way, the full propagator is constructed to describe exclusively the diffusive dynamics in the bulk. As the stochastic solution of the Skorokhod equation (1) is unique [10, 11], the full propagator is well-defined.

In order to incorporate surface reactions, we exploit the microscopic interpretation of the boundary condition (8) by introducing a thin reactive layer of width $a$ near the boundary and counting the number of times, $\mathcal{N}_{t}^{a}$, that the diffusing particle has crossed this layer up to time $t$. Each crossing can be understood as an encounter with the boundary during which the particle might interact with it. The self-similar character of Brownian motion implies that after the first encounter with the boundary, the particle returns infinitely many times to that boundary [42]. As a consequence, the number $\mathcal{N}_{t}^{a}$ diverges as $a\to 0$, but its rescaling by $a$ yields a nontrivial limit – the boundary local time:

This representation is equivalent to Eq. (2) and yields an alternative approximation to the boundary local time.

At each encounter, the particle can react with the boundary with some probability $p_{a}\simeq a\kappa/D$ or resume its diffusion in the bulk [31, 33, 43]. Such reaction attempts are supposed to be independent from each other and from the diffusion process. Denoting by ${\mathcal{T}}$ the time of the successful reaction event, one gets in the limit $a\to 0$:

where we introduced an independent random threshold $\hat{\ell}$ obeying the exponential law:

In other words, the successful reaction event occurs when the boundary local time $\ell_{t}$ crosses the random threshold $\hat{\ell}$. As first suggested in [35], this relation can actually be used as the definition of the first-reaction time ${\mathcal{T}}$:

in analogy with the definition of the first-passage time: ${\mathcal{T}}_{\rm FPT}=\inf\{t>0~{}:~{}\bm{X}_{t}\in\partial\Omega\}$. Moreover, as the boundary local time remains zero until the first encounter, the first-passage time can also be defined as ${\mathcal{T}}_{\rm FPT}=\inf\{t>0~{}:~{}\ell_{t}>0\}$, i.e., as the first crossing of the zero threshold. The boundary local time offers thus a unified way of treating perfectly and partially reactive boundaries.

By definition, the conventional propagator $G_{q}(\bm{x},t|\bm{x}_{0})$ is the probability density of finding the particle in a vicinity of $\bm{x}$ at time $t$, given that the particle has started at $\bm{x}_{0}$ and not reacted on the boundary up to time $t$. According to Eq. (13), the survival of the particle means that $\ell_{t}<\hat{\ell}$. As a consequence, the conventional propagator is given by the full propagator $P(\bm{x},\ell,t|\bm{x}_{0})$, multiplied by the survival probability of that particle up to time $t$, ${\mathbb{P}}\{{\mathcal{T}}>t\}={\mathbb{P}}\{\ell_{t}<\hat{\ell}\}$, and integrated over all possible values of the boundary local time:

The parameter $q$ stands in the subscript of the propagator to highlight its dependence on $q$ through the boundary condition (8). In other words, the propagator for partially reactive boundary is obtained as the Laplace transform (with respect to $\ell$) of the full propagator for purely reflecting boundary. This relation was established in [8] for ordinary diffusion without drift. Whenever surface reactions are independent from bulk dynamics, the above arguments can be applied, in particular, Eq. (16) holds for more general diffusion processes with a drift. Note also that this relation can also be deduced from the rigorous probabilistic analysis of a more general Robin boundary value problem presented in [44]. We emphasize that the boundary local time $\ell_{t}$ depends exclusively on the diffusive dynamics in the confining domain $\Omega$ with a purely reflecting boundary, whereas the threshold $\hat{\ell}$ depends exclusively on the reactivity $q$, allowing one to disentangle these two aspects of diffusion-controlled reactions. This disentanglement opens a way to investigate much more elaborate surface reaction mechanisms when the exponential distribution of the random threshold $\hat{\ell}$ is replaced by another distribution (see [8] for details). As a consequence, the unique full propagator $P(\bm{x},\ell,t|\bm{x}_{0})$ can determine a variety of conventional propagators by choosing an appropriate stopping condition. In summary, the stochastic foundation (1) for both conventional and encounter-based approaches is the same, but the ways of incorporating surface reactions are different (Robin boundary condition versus random threshold $\hat{\ell}$).

While the full propagator $P(\bm{x},\ell,t|\bm{x}_{0})$ determines any conventional propagator $G_{q}(\bm{x},t|\bm{x}_{0})$ via Eq. (16), its inverse Laplace transform with respect to $q$ can in turn be used to access the full propagator. However, the implicit dependence of the conventional propagator on $q$ as the parameter in the boundary condition (8) makes difficult further explorations of such an inverse, even numerically. To overcome this limitation, a spectral representation of the full propagator was derived in [8] for ordinary diffusion. Here we aim at extending this representation to restricted diffusion with a gradient drift.

In the following, we mainly deal with Laplace-transformed quantities with respect to time $t$, denoted by tilde. For instance,

is the Green’s function for the Fokker-Planck equation, which admits both forward and backward forms:

and

where

is the Laplace-transformed probability flux density.

In [8], we considered ordinary restricted diffusion (without drift) and employed the so-called Dirichlet-to-Neumann operator ${\mathcal{M}}_{p}$ that associates to a given function $f$ on the boundary $\partial\Omega$ another function on that boundary such that ${\mathcal{M}}_{p}f=(\partial_{n}u)_{|\partial\Omega}$, where $u(\bm{x})$ satisfies $(p-D\Delta)u(\bm{x})=0$ in $\Omega$ and $u_{|\partial\Omega}=f$. In other words, ${\mathcal{M}}_{p}$ maps the Dirichlet boundary condition onto Neumann boundary condition. Relying on the fact that the self-adjoint operator ${\mathcal{M}}_{p}$ on a bounded boundary has a discrete spectrum (see [45, 46, 47, 48, 49, 50, 51] for mathematical details), we derived the following spectral decomposition

where $\mu_{n}^{(p)}$ are the eigenvalues of ${\mathcal{M}}_{p}$, while $V_{n}^{(p)}(\bm{x})$ are projections of $\tilde{j}_{\infty}(\bm{x},p|\bm{x}_{0})$ onto the eigenfunctions $v_{n}^{(p)}$ of that operator (see below for details; note that Eq. (21) is not used in the following derivation and is reproduced here just for motivation). The inverse Laplace transform of this relation with respect to $q$ yielded the spectral decomposition of the Laplace-transformed full propagator:

A naive extension of this approach to restricted diffusion with drift would fail. In fact, even though an appropriate Dirichlet-to-Neumann operator could be introduced by replacing the equation $(p-D\Delta)u=0$ by $(p-{\mathcal{L}}_{\bm{x}})u=0$ or $(p-{\mathcal{L}}^{\dagger}_{\bm{x}_{0}})u=0$, such an extension would not in general be self-adjoint and thus would require much more elaborate spectral tools.

To overcome this limitation, we reformulate the problem in terms of the symmetrized Fokker-Planck operator [2]

where the second term in the last relation is the multiplication by the explicit drift-dependent function (here, the operator $\nabla$ acts only on $\bm{\mu}(\bm{x})$). To avoid mathematical subtleties, we assume that the second term is a smooth bounded function on $\Omega$. With the aid of this expression, the Fokker-Planck operator and its adjoint can be written as

In order to reduce the original problem to that with $\bar{{\mathcal{L}}}$, one can represent the Green’s function as

so that forward and backward equations imply two equivalent boundary value problems for the new function $\tilde{g}_{q}(\bm{x},p|\bm{x}_{0})$:

and

with

As the forward and backward problems are now identical, their solution is symmetric with respect to the exchange of $\bm{x}$ and $\bm{x}_{0}$: $\tilde{g}_{q}(\bm{x},p|\bm{x}_{0})=\tilde{g}_{q}(\bm{x}_{0},p|\bm{x})$.

We search for the following representation:

where the new function $u(\bm{x},p|\bm{x}_{0})$ satisfies

subject to the boundary condition at $\bm{x}_{0}\in\partial\Omega$:

i.e.,

where

Since the normal derivative $\partial_{n_{0}}$ acts on $\bm{x}_{0}$ of $\tilde{G}_{q}(\bm{x},p|\bm{x}_{0})$, this function differs in general from $\tilde{j}_{\infty}(\bm{x},p|\bm{x}_{0})$ given in Eq. (20). We also emphasize that the order of $\bm{x}_{0}$ and $\bm{x}$ in this notation has been changed, in order to keep the convention that the first argument of $j$ is a point on the boundary.

In the next step, we introduce a Dirichlet-to-Neumann operator ${\mathcal{M}}_{p}$ that associates to a given function $f$ on $\partial\Omega$ another function on $\partial\Omega$ such that

The operator $\bar{{\mathcal{L}}}$ in Eq. (23) having the form of a Schrödinger operator, is self-adjoint; in addition, as the domain is bounded, the drift-dependent term is a bounded perturbation to the Laplace operator. In analogy with the Laplacian case (see [47, 48] for mathematical details), one can expect that, under mild assumptions on the drift term, the Dirichlet-to-Neumann operator ${\mathcal{M}}_{p}$ is self-adjoint, with a discrete spectrum of real eigenvalues $\{\mu_{n}^{(p)}\}$, while its eigenfunctions $\{v_{n}^{(p)}\}$ form a complete orthonormal basis in $L_{2}(\partial\Omega)$. As a proof of this mathematical statement is beyond the scope of the present paper, we will use it as a conjecture in the following presentation. In other words, we focus on such drifts $\bm{\mu}(\bm{x})$, for which this statement is valid. The Dirichlet-to-Neumann operator ${\mathcal{M}}_{p}$ allows us to rewrite the boundary condition (34) in an operator form, from which

Finally, one needs to extend $u(\bm{x},p|\bm{x}_{0})$ to any $\bm{x}_{0}\in\Omega$. For this purpose, one multiplies Eq. (27) with $q=\infty$ by $u(\bm{x},p|\bm{x}^{\prime})$ and subtracts from it the equation $(p-\bar{{\mathcal{L}}}_{\bm{x}})u(\bm{x},p|\bm{x}^{\prime})=0$ multiplied by $\tilde{g}_{\infty}(\bm{x}^{\prime},p|\bm{x}_{0})$. The integration of both parts over $\bm{x}^{\prime}$ yields

where we used Green’s relations and the Dirichlet boundary condition $\tilde{G}_{\infty}(\bm{x},p|\bm{x}_{0})=0$ at $\bm{x}\in\partial\Omega$. As a consequence, we get

Using the eigenfunctions of ${\mathcal{M}}_{p}$, one gets the spectral decomposition of the Green’s function:

where

When there is no drift, one has $\Phi(\bm{s})=0$ and thus $V^{{}^{\prime}(p)}_{n}(\bm{x})=V^{(p)}_{n}(\bm{x})$, so that Eq. (39) is reduced to the former result (21) derived in [8]. The Laplace inversion of Eq. (39) with respect to $q$ yields

This is our main theoretical result, which generalizes the former relation (22) to the case of restricted diffusion with a gradient drift.

As said earlier, the full propagator determines various characteristics of diffusion-influenced reactions. For instance, integrating Eqs. (39, 42) over $\bm{x}\in\Omega$, one accesses the Laplace-transformed survival probability and the Laplace-transformed probability density of the boundary local time, respectively:

and

where $\circ$ denotes the marginalized variable $\bm{x}$, and

In A, we derive the following representation for this quantity:

These expressions generalize our former results for ordinary diffusion without drift from [52, 53].

We recall that the survival probability determines the distribution of the first-reaction time ${\mathcal{T}}$ on a partially reactive boundary [4]: $S_{q}(t|\bm{x}_{0})={\mathbb{P}}_{\bm{x}_{0}}\{{\mathcal{T}}>t\}$. At the same time, $S_{q}(t|\bm{x}_{0})={\mathbb{E}}_{\bm{x}_{0}}\{e^{-q\ell_{t}}\}$ is the generating function of the boundary local time $\ell_{t}$ [8, 41] that allows one to compute its moments as

from which

Multiplying both sides by $p$, one can also interpret the integral over $t$ as the double average of the random variable $\ell_{\tau}^{k}$ with the random exponentially distributed stopping time $\tau$. Similarly, $p\tilde{P}(\circ,\ell,p|\bm{x}_{0})$ can be interpreted as the probability density of $\ell_{\tau}$:

Setting $q=0$ in Eq. (43) and noting that $\tilde{S}_{0}(p|\bm{x}_{0})=1/p$ (since $S_{0}(t|\bm{x}_{0})=1$ for a fully inert boundary), one gets

Substituting this expression back into Eq. (43), one deduces an equivalent spectral decomposition:

from which one also derives the Laplace-transformed probability density of the first-reaction time:

As $H_{q}(t|\bm{x}_{0})$ can be interpreted as the probability flux of particles, started from $\bm{x}_{0}$, onto the reactive boundary, the total flux $J_{q}(t)$ is obtained by averaging $H_{q}(t|\bm{x}_{0})$ with the initial concentration $c(\bm{x}_{0})$ of these particles. If the initial concentration is uniform, $c(\bm{x}_{0})=c_{0}$, one gets in the Laplace domain: