# An encounter-based approach to the escape problem

Denis S. Grebenkov<EMAIL_ADDRESS>Laboratoire de Physique
de la Matière Condensée,
CNRS – Ecole Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau,
France

###### Abstract

We revise the encounter-based approach to imperfect diffusion-controlled
reactions, which employs the statistics of encounters between a diffusing
particle and the reactive region to implement surface reactions. We extend
this approach to deal with a more general setting, in which the reactive
region is surrounded by a reflecting boundary with an escape region. We derive
a spectral expansion for the full propagator and investigate the behavior and
probabilistic interpretations of the associated probability flux density. In
particular, we obtain the joint probability density of the escape time and the
number of encounters with the reactive region before escape, and the
probability density of the first-crossing time of a prescribed number of
encounters. We briefly discuss generalizations of the conventional Poissonian-
type surface reaction mechanism described by Robin boundary condition and
potential applications of this formalism in chemistry and biophysics.

Diffusion-controlled reaction, escape problem, first-passage time, surface
reaction, boundary local time, encounters, Robin boundary condition, reflected
Brownian motion, restricted diffusion

###### pacs:

02.50.-r, 05.40.-a, 02.70.Rr, 05.10.Gg

##  I Introduction

Diffusion-controlled reactions play an important role for various chemical and
biophysical applications Alberts ; Rice ; Lauffenburger ; Schuss ; Lindenberg
; Bressloff13 . In a typical setting, a particle diffuses inside a confining
domain $\Omega$ toward a target region, on which it can react or trigger a
specific event. One can think of a protein searching for a specific site on a
DNA molecule, or a molecule in a chemical reactor searching for a catalytic
germ to be transformed. Such diffusion-controlled reactions are often
described in terms of the first-passage time to the target or, more generally,
to the reaction event Redner ; Metzler ; Condamin07 ; Benichou10 ; Holcman13 ;
Benichou14 ; Grebenkov16 ; Guerin16 ; Lanoiselee18 . The diffusive dynamics of
a single molecule is usually characterized by a propagator,
$G_{q}(\bm{x},t|\bm{x}_{0})$, i.e., the probability density that a molecule
started from $\bm{x}_{0}$ at time $0$ has arrived in a vicinity of point
$\bm{x}$ at a later time $t$, without being reacted Redner ; Gardiner ;
VanKampen . For ordinary diffusion, the propagator satisfies the diffusion (or
heat) equation,

$\partial_{t}G_{q}(\bm{x},t|\bm{x}_{0})=D\Delta
G_{q}(\bm{x},t|\bm{x}_{0})\quad(\bm{x}\in\Omega),$ (1)

subject to the initial condition
$G_{q}(\bm{x},0|\bm{x}_{0})=\delta(\bm{x}-\bm{x}_{0})$ with a Dirac
distribution $\delta(\bm{x}-\bm{x}_{0})$, where $D$ is the diffusion
coefficient of the molecule, and $\Delta$ is the Laplace operator (the meaning
of the subscript $q$ will be explained below). In turn, the reactive
properties of the boundary of the confining domain are usually incorporated
through the Robin boundary condition,

$-D\partial_{n}G_{q}(\bm{x},t|\bm{x}_{0})=\kappa(\bm{x})\,G_{q}(\bm{x},t|\bm{x}_{0})\quad(\bm{x}\in\partial\Omega),$
(2)

where $\partial_{n}$ is the normal derivative on the boundary, oriented
outwards the confining domain, and $\kappa(\bm{x})$ is the reactivity at a
boundary point $\bm{x}$. This condition, which was put forward in chemical
physics by Collins and Kimball Collins49  and broadly employed afterwards
Berg77 ; Sano79 ; Brownstein79 ; Weiss86 ; Powles92 ; Sapoval94 ; Sapoval02 ;
Grebenkov05 ; Traytak07 ; Singer08 ; Bressloff08 ; Grebenkov10 ; Lawley15 ;
Grebenkov15 ; Serov16 ; Bressloff17 ; Grebenkov17 ; Piazza19 , postulates that
the diffusive flux from the bulk (the left-hand side) is proportional to the
propagator on the boundary. Depending on the type of surface reaction, the
proportionality coefficient $\kappa(\bm{x})$ (in units m/s) is called
reactivity, permeability, or surface relaxivity. In the context of bimolecular
reactions, it can also be related to the forward reaction constant. The
reactivity $\kappa(\bm{x})$ can range from $0$ for inert impermeable boundary
to $+\infty$ for a perfectly reactive boundary, on which the reaction occurs
upon the first arrival of the particle onto that boundary.

A space-dependent reactivity $\kappa(\bm{x})$ allows one to implement
heterogeneous reactivity patterns on a catalytic surface or to describe in a
unified way the effects of reactive targets and restricting inert boundaries
of a porous medium. However, theoretical description of such general
diffusion-controlled reactions is rather challenging (see Grebenkov19b  and
references therein). For this reason, one often focuses on a simpler yet
relevant setting of a constant or piecewise constant reactivity. If
$\kappa(\bm{x})$ is constant, one can employ standard spectral expansions over
the eigenfunctions of the Laplace operator Redner ; Gardiner ; VanKampen .
Moreover, for simple confining domains such as spheres or parallelepipeds,
these eigenfunctions are known exactly that facilitates the analysis of
diffusion-controlled reactions Carslaw ; Crank ; Thambynayagam ; Grebenkov13 .
Another common setting is the case of a reactive target surrounded by a
reflecting boundary so that $\kappa(\bm{x})=\kappa$ on the target surface, and
$\kappa(\bm{x})=0$ on the reflecting boundary. In this case, one deals with
mixed Robin-Neumann (or Dirichlet-Neumann for $\kappa=\infty$) boundary
conditions. When the target is small, one can employ powerful asymptotic tools
to approximate various quantities such as, for instance, the mean reaction
time on the target or the decay rate of the survival probability Mazya85 ;
Ward93 ; Kolokolnikov05 ; Singer06a ; Singer06b ; Singer06c ; Schuss07 ;
Benichou08 ; Pillay10 ; Cheviakov10 ; Cheviakov11 ; Cheviakov12 ; Holcman14 ;
Agranov18 ; Grebenkov19c . Moreover, if the target is interpreted as a “hole”
in the otherwise reflecting impenetrable boundary, one speaks about the escape
or exit problem. As the reaction is understood here as an escape event, the
reaction time is called the escape or exit time. Finally, if there are many
targets, their competition for the diffusing particle can be characterized by
splitting probabilities Traytak92 ; Condamin06 ; Chevalier11 ; Galanti16 ;
Grebenkov19f ; Grebenkov20f ; Klinger22 .

In this paper, we consider a more general situation, in which the boundary
$\partial\Omega$ of the confining domain is split into three disjoint parts,

$\partial\Omega=\partial\Omega_{R}\cup\partial\Omega_{N}\cup\partial\Omega_{D},$
(3)

which represent the reactive target $\partial\Omega_{R}$, the inert reflecting
boundary $\partial\Omega_{N}$ and the escape region $\partial\Omega_{D}$ (Fig.
1). In this way, one can describe an important class of diffusion-controlled
reactions, in which the diffusing particle can leave the confining domain
through an escape region or be destroyed on it, without being reacted on the
target region. This is a common setting for many biochemical reactions inside
a living cell; for instance, proteins can be disassembled before finding their
receptors, while ions can leave the cytoplasm through the ion channels on the
plasma membrane. We assume that the particle disappears immediately after the
first arrival on the escape region $\partial\Omega_{D}$. In mathematical
terms, such a composed boundary can be implemented through the mixed Robin-
Dirichlet-Neumann boundary conditions:

$\displaystyle\partial_{n}G_{q}(\bm{x},t|\bm{x}_{0})+qG_{q}(\bm{x},t|\bm{x}_{0})$
$\displaystyle=0\quad(\bm{x}\in\partial\Omega_{R}),$ (4a) $\displaystyle
G_{q}(\bm{x},t|\bm{x}_{0})$
$\displaystyle=0\quad(\bm{x}\in\partial\Omega_{D}),$ (4b)
$\displaystyle\partial_{n}G_{q}(\bm{x},t|\bm{x}_{0})$
$\displaystyle=0\quad(\bm{x}\in\partial\Omega_{N}),$ (4c)

where $q=\kappa/D$ is proportional to the (constant) reactivity $\kappa$ of
the target region $\partial\Omega_{R}$ (note that subscripts of
$\partial\Omega_{R}$, $\partial\Omega_{N}$ and $\partial\Omega_{D}$ refer to
the corresponding type of boundary condition). Alternatively, these conditions
can describe a two-target problem: a partially reactive target
$\partial\Omega_{R}$ and a perfectly reactive target $\partial\Omega_{D}$
(“reaction” on $\partial\Omega_{D}$ is interpreted here as the escape event).
However, as we focus on reactions on the target region $\partial\Omega_{R}$,
we keep speaking about a single-target problem in the presence of escape
events.

Figure 1:  Illustration of a confining domain $\Omega$ whose boundary
$\partial\Omega$ is split into three disjoint parts: a target region
$\partial\Omega_{R}$ (the boundary of a gray obstacle, e.g., a catalytic
germ), a reflecting region $\partial\Omega_{N}$ (in green), and an escape
region $\partial\Omega_{D}$ (in red). A simulated trajectory of a particle
that started from a point $\bm{x}_{0}$ and diffused until its escape, is shown
in magenta. This particle is always reflected from $\partial\Omega_{N}$ but
may either be reflected from or react on $\partial\Omega_{R}$.

In order to solve this problem, one could still introduce the Laplacian
eigenfunctions that satisfy the same mixed boundary conditions. Here we follow
an alternative way and generalize the encounter-based approach that was
developed in Grebenkov20  for the particular case of a constant reactivity on
the whole boundary (when $\partial\Omega=\partial\Omega_{R}$ and
$\partial\Omega_{N}=\partial\Omega_{D}=\emptyset$). This approach relies on
the concept of the boundary local time $\ell_{t}$, which can represent a
rescaled number of encounters between the diffusing particle and the boundary
up to time $t$. In this way, one can first investigate the diffusive dynamics
inside a confining domain with reflecting boundary and then implement surface
reactions explicitly (see below). Moreover, one can go beyond the conventional
Poissonian-type reactions described by the Robin boundary condition (2) and
implement saturation or activation effects or, more generally, encounter-
dependent reactivity Grebenkov20 . Our generalization allows one to
investigate the effects of the escape event onto such reactions. We also
provide probabilistic interpretations of the probability flux density that
were not discussed enough in earlier works. In particular, we obtain the joint
probability density of the position, boundary local time and the escape time.

The paper is organized as follows. In Sec. II, we present the main theoretical
results. After introducing the necessary elements of the conventional and
encounter-based approaches in Sec. II.1 and Sec. II.2, we discuss the
distribution of the boundary local time (Sec. II.3) and restrictions of the
probability flux density to $\partial\Omega_{D}$ and $\partial\Omega_{R}$
(Secs. II.4 and II.5). In Sec. III, we illustrate the behavior of the derived
quantities for a particle diffusing between two concentric spheres. In this
case, all the “ingredients” of the encounter-based approach can be found
explicitly. Section IV presents further discussions, conclusions and
perspectives.

##  II Main results

###  II.1 Conventional approach

In order to highlight the advantages of the encounter-based approach, we
briefly recall several characteristics of diffusion-controlled reactions
accessed within the conventional approach that relies on the propagator
$G_{q}(\bm{x},t|\bm{x}_{0})$. We consider a point-like particle diffusing in a
bounded confining Euclidean domain $\Omega\subset\mathbb{R}^{d}$ with smooth
boundary $\partial\Omega$, which is partitioned into reactive
($\partial\Omega_{R}$), reflecting ($\partial\Omega_{N}$) and escape
($\partial\Omega_{D}$) parts. As the diffusing particle can disappear due to
either escape through $\partial\Omega_{D}$ or reaction on
$\partial\Omega_{R}$, one can naturally introduce two first-passage times: the
escape time

$T_{q}=\inf\\{t>0~{}:~{}\bm{X}_{t}\in\partial\Omega_{D}\\}$ (5)

as the first-passage time to the escape region $\partial\Omega_{D}$, and the
reaction time $\tau_{q}$ as the first instance of the reaction on
$\partial\Omega_{R}$ (its formal definition is nontrivial and will be given in
Sec. II.5). The subscript $q$ highlights that the distributions of both random
variables $T_{q}$ and $\tau_{q}$ depend on the reactivity of the target region
because they are determined by the probability flux density,

$j_{q}(\bm{x},t|\bm{x}_{0})=-D\partial_{n}G_{q}(\bm{x},t|\bm{x}_{0})\quad(\bm{x}\in\partial\Omega).$
(6)

In fact, the restriction of $j_{q}(\bm{x},t|\bm{x}_{0})$ to the escape region
$\partial\Omega_{D}$ is the joint probability density of the escape location
$\bm{X}_{T_{q}}$ and its time $T_{q}$. In turn, the restriction of
$j_{q}(\bm{x},t|\bm{x}_{0})$ to the target region $\partial\Omega_{R}$ is the
joint probability density of the reaction location $\bm{X}_{\tau_{q}}$ and its
time $\tau_{q}$ (note that the restriction of $j_{q}(\bm{x},t|\bm{x}_{0})$ to
$\partial\Omega_{N}$ is strictly zero due to the Neumann boundary condition
(4c)). If the position does not matter, it can be averaged out to get the
(marginal) probability densities of the escape time $T_{q}$ and of the
reaction time $\tau_{q}$:

$\displaystyle J_{q}^{D}(t|\bm{x}_{0})$
$\displaystyle=\int\limits_{\partial\Omega_{D}}d\bm{x}\,j_{q}(\bm{x},t|\bm{x}_{0}),$
(7a) $\displaystyle J_{q}^{R}(t|\bm{x}_{0})$
$\displaystyle=\int\limits_{\partial\Omega_{R}}d\bm{x}\,j_{q}(\bm{x},t|\bm{x}_{0}).$
(7b)

As the particle can either escape or react, none of these densities is
normalized to $1$; in turn, one has

$\int\limits_{0}^{\infty}dt\,\bigl{[}J_{q}^{D}(t|\bm{x}_{0})+J_{q}^{R}(t|\bm{x}_{0})\bigr{]}=1.$
(8)

To show this normalization, one can integrate the diffusion equation (1) over
$\bm{x}\in\Omega$ and use the Green’s formula and mixed boundary conditions
(4) to get the continuity equation

$\partial_{t}S_{q}(t|\bm{x}_{0})=-\bigl{[}J_{q}^{D}(t|\bm{x}_{0})+J_{q}^{R}(t|\bm{x}_{0})\bigr{]},$
(9)

where

$S_{q}(t|\bm{x}_{0})=\int\limits_{\Omega}d\bm{x}\,G_{q}(\bm{x},t|\bm{x}_{0})$
(10)

is the survival probability of the particle up to time $t$. The integral of
Eq. (9) over $t$ from $0$ to infinity yields the normalization (8), given that
$S_{q}(0|\bm{x}_{0})=1$ and $S_{q}(t|\bm{x}_{0})\to 0$ as $t\to\infty$ for
diffusion in a bounded domain.

###  II.2 Encounter-based approach

In the encounter-based approach, one first characterizes purely diffusive
dynamics inside a bounded confining domain $\Omega$ with a reflecting boundary
$\partial\Omega$, and then incorporates surface reactions on $\partial\Omega$,
i.e., transforms the reflecting boundary into the reactive one Grebenkov20 .
For this purpose, one uses the boundary local time $\ell_{t}$, which was first
introduced by Lévy Levy  and then extensively employed in mathematical
literature on stochastic processes Ito ; Freidlin . On one hand, the boundary
local time can be defined as the renormalized residence time near the boundary
$\partial\Omega$:

$\ell_{t}=\lim\limits_{a\to
0}\frac{D}{a}\int\limits_{0}^{t}dt^{\prime}\,\Theta(a-|\bm{X}_{t^{\prime}}-\partial\Omega|),$
(11)

where $|\bm{X}_{t^{\prime}}-\partial\Omega|$ is the Euclidean distance between
the position $\bm{X}_{t^{\prime}}$ of the particle at time $t^{\prime}$ and
the boundary $\partial\Omega$, and $\Theta(z)$ is the Heaviside function:
$\Theta(z)=1$ for $z>0$ and $0$ otherwise. The integral in Eq. (11) defines
the residence (also known as occupation or sojourn) time in a thin layer of
width $a$ near the boundary. As the boundary $\partial\Omega$ has a lower
dimension as compared to the domain $\Omega$, the residence time vanishes as
the layer shrinks (when $a\to 0$); in turn, its renormalization by $a$ yields
the nontrivial limit (11). The boundary local time should not be confused with
the local time in a bulk point, which has been intensively studied, especially
for diffusive processes in one dimension Borodin ; Majumdar05 . Despite its
name, the boundary local time has units of length. On the other hand, the
boundary local time can be equivalently defined as a rescaled limit of the
number $\mathcal{N}_{t}^{a}$ of downcrossings of the boundary layer up to time
$t$: $\ell_{t}=\lim\limits_{a\to 0}a\mathcal{N}_{t}^{a}$. As each downcrossing
can be interpreted as an encounter of the particle with the boundary, the
boundary local time $\ell_{t}$ characterizes the statistics of such encounters
Grebenkov21 . The diffusive dynamics can then be described either by a
stochastic differential equation for the pair $(\bm{X}_{t},\ell_{t})$, or by
the so-called full propagator $P(\bm{x},\ell,t|\bm{x}_{0})$, i.e., the joint
probability density of getting the values $(\bm{x},\ell)$ for the pair
$(\bm{X}_{t},\ell_{t})$. Once the full propagator is known, one can implement
various surface reaction mechanisms (see Sec. II.5).

By construction, the boundary local time $\ell_{t}$ characterizes encounters
with the whole boundary $\partial\Omega$ that does not allow one to implement
different reaction mechanisms on different subsets of the boundary. This
problem has been discussed and partly resolved in Grebenkov20c  by introducing
a proper boundary local time $\ell_{t}^{i}$ on each subset of interest. In
this case, one would deal with the joint probability density for $\bm{X}_{t}$
and boundary local times on all these subsets. Even though a formal way for
computing this density was proposed in Grebenkov20c , its practical
implementation was realized only for simple geometric settings (e.g., an
interval).

In the context of the escape problem that we consider in this paper, one can
follow a different strategy. As we are interested in describing reactions
exclusively on the target region $\partial\Omega_{R}$, we need to know the
statistics of encounters with that particular region. In other words, one can
modify the above definition of the boundary local time to count encounters
only with $\partial\Omega_{R}$:

$\ell_{t}^{R}=\lim\limits_{a\to
0}\frac{D}{a}\int\limits_{0}^{t}dt^{\prime}\,\Theta(a-|\bm{X}_{t^{\prime}}-\partial\Omega_{R}|).$
(12)

This relation defines a non-decreasing stochastic process $\ell_{t}^{R}$
(starting from $\ell_{0}^{R}=0$) that increases at each encounter with
$\partial\Omega_{R}$. In the same vein, we introduce the full propagator
$P(\bm{x},\ell,t|\bm{x}_{0})$ as the joint probability density for
$\bm{X}_{t}$ and $\ell_{t}^{R}$ (not $\ell_{t}$), under the condition of no
escape through $\partial\Omega_{D}$ up to time $t$. As previously, the target
region $\partial\Omega_{R}$ is treated at this step as reflecting, i.e., the
particle described by the full propagator can disappear only on the escape
region $\partial\Omega_{D}$.

At the next step, one can introduce reaction events on $\partial\Omega_{R}$
following the probabilistic arguments from Grebenkov20 . For this purpose, one
can consider a thin reactive layer of width $a$ near the target region
$\partial\Omega_{R}$. Once the particle enters this layer, surface reaction
may be described by a standard first-order reaction kinetics, with the rate
$k=\kappa/a$. Since the residence time of the particle within this layer up to
time $t$ is approximately equal to $\ell a/D$ according to Eq. (12), the
probability of no reaction on $\partial\Omega_{R}$ up to $t$ is then
$e^{-k(\ell a/D)}=e^{-q\ell}$, where $q=\kappa/D$. As a consequence, one
deduces the following relation between the conventional and full propagators:

$G_{q}(\bm{x},t|\bm{x}_{0})=\int\limits_{0}^{\infty}d\ell\,e^{-q\ell}\,P(\bm{x},\ell,t|\bm{x}_{0}).$
(13)

On the left-hand side, the conventional propagator
$G_{q}(\bm{x},t|\bm{x}_{0})$ satisfying Eqs. (1, 4), describes diffusion from
$\bm{x}_{0}$ to $\bm{x}$ in time $t$, under the condition of no escape through
$\partial\Omega_{D}$ and no reaction on $\partial\Omega_{R}$ up to $t$. On the
right-hand side, $P(\bm{x},\ell,t|\bm{x}_{0})$ describes diffusion from
$\bm{x}_{0}$ to $\bm{x}$ in time $t$, under the condition of no escape through
$\partial\Omega_{D}$ and of getting the boundary local time $\ell_{t}^{R}$
equal to $\ell$. In turn, the factor $e^{-q\ell}$ incorporates the probability
of no reaction on $\partial\Omega_{R}$ for any realized value $\ell$ of the
boundary local time $\ell_{t}^{R}$, while the integral over $\ell$ sums up
contributions from all possible realizations.

The relation (13) plays the central role in this work. On one hand, it allows
one to incorporate the conventional surface reactions as the Laplace transform
of the full propagator $P(\bm{x},\ell,t|\bm{x}_{0})$ with respect to $\ell$.
Importantly, the reactivity parameter $q$ enters explicitly through the
exponential factor $e^{-q\ell}$, whereas it appeared implicitly in the
conventional description as a coefficient in the Robin boundary condition
(4a). Finally, one can replace the factor $e^{-q\ell}$, which is reminiscent
for an exponential probability law, by another function, allowing one to
implement various surface reaction mechanisms (see Sec. II.5). On the other
hand, the inverse Laplace transform of the conventional propagator
$G_{q}(\bm{x},t|\bm{x}_{0})$ with respect to $q$ gives access to the full
propagator:

$P(\bm{x},\ell,t|\bm{x}_{0})=\mathcal{L}_{q\to\ell}^{-1}\bigl{\\{}G_{q}(\bm{x},t|\bm{x}_{0})\bigr{\\}}.$
(14)

Unfortunately, an implicit dependence of $G_{q}(\bm{x},t|\bm{x}_{0})$ on $q$
often prevents using this inversion and thus urges for another representation
for the full propagator.

To achieve this goal, we extend the spectral expansion of the full propagator
developed in Grebenkov20 . For this purpose, we introduce an extension of the
so-called Dirichlet-to-Neumann operator $\mathcal{M}_{p}$ that associates to a
given function $f$ on the target region $\partial\Omega_{R}$ another function
$g$ on that region such that

$\mathcal{M}_{p}f=g=\left.(\partial_{n}u)\right|_{\partial\Omega_{R}},$ (15)

where $u$ is the solution of the following boundary value problem with $p\geq
0$:

$\displaystyle(p-D\Delta)u$ $\displaystyle=0\quad(\bm{x}\in\Omega),$ (16a)
$\displaystyle u$ $\displaystyle=f\quad(\bm{x}\in\partial\Omega_{R}),$ (16b)
$\displaystyle u$ $\displaystyle=0\quad(\bm{x}\in\partial\Omega_{D}),$ (16c)
$\displaystyle\partial_{n}u$
$\displaystyle=0\quad(\bm{x}\in\partial\Omega_{N}).$ (16d)

In other words, the operator $\mathcal{M}_{p}$ transforms the Dirichlet
boundary condition $u=f$ on $\partial\Omega_{R}$ into an equivalent Neumann
boundary condition $\partial_{n}u=g$ on $\partial\Omega_{R}$, keeping
unchanged the Dirichlet and Neumann conditions on $\partial\Omega_{D}$ and
$\partial\Omega_{N}$ respectively. While the conventional Dirichlet-to-Neumann
operator acted on functions on the whole boundary $\partial\Omega$, our
extension acts on functions on the subset $\partial\Omega_{R}$ of the
boundary. The spectral properties of the conventional Dirichlet-to-Neumann
operator have been intensively studied in mathematical literature Arendt14 ;
Daners14 ; terElst14 ; Behrndt15 ; Arendt15 ; Hassell17 ; Girouard17 . Most of
these properties are expected to be valid for our extension so that
$\mathcal{M}_{p}$ is a pseudo-differential self-adjoint operator. Since the
target region $\partial\Omega_{R}$ is bounded, the spectrum of
$\mathcal{M}_{p}$ is discrete, with an infinite set of positive eigenvalues
$\mu_{k}^{(p)}$ ($k=0,1,2,\ldots$), that can be enumerated in the increasing
order:

$0\leq\mu_{0}^{(p)}\leq\mu_{1}^{(p)}\leq\ldots\leq\mu_{k}^{(p)}\leq\ldots\nearrow+\infty.$
(17)

In turn, the associated eigenfunctions $v_{k}^{(p)}(\bm{x})$ form a complete
orthonormal basis of the functional space $L_{2}(\partial\Omega_{R})$ of
square integrable functions on $\partial\Omega_{R}$. A rigorous formulation
and demonstration of these mathematical properties are beyond the scope of
this paper. From the mathematical point of view, one can consider them as
conjectural extensions of the well-known conventional case.

The eigenbasis of the Dirichlet-to-Neumann operator can serve for getting a
spectral expansion of the full propagator. Skipping technical details (given
in Grebenkov20 ), we sketch here the main steps of this derivation. In the
first step, the Laplace transform of the conventional propagator
$G_{q}(\bm{x},t|\bm{x}_{0})$ with respect to time $t$,

$\tilde{G}_{q}(\bm{x},p|\bm{x}_{0})=\int\limits_{0}^{\infty}dt\,e^{-pt}\,G_{q}(\bm{x},t|\bm{x}_{0}),$
(18)

reduces the diffusion equation (1) to the inhomogeneous modified Helmholtz
equation:

$(p-D\Delta)\tilde{G}_{q}(\bm{x},p|\bm{x}_{0})=\delta(\bm{x}-\bm{x}_{0}),$
(19)

subject to the same boundary condition (here and below, tilde denotes Laplace-
transformed quantities with respect to $t$). Writing
$\tilde{G}_{q}(\bm{x},p|\bm{x}_{0})=\tilde{G}_{\infty}(\bm{x},p|\bm{x}_{0})+\tilde{g}_{q}(\bm{x},p|\bm{x}_{0})$
with an unknown function $\tilde{g}_{q}(\bm{x},p|\bm{x}_{0})$, one eliminates
$\delta(\bm{x}-\bm{x}_{0})$ from the right-hand side. As
$\tilde{g}_{q}(\bm{x},p|\bm{x}_{0})$ satisfies
$(p-D\Delta)\tilde{g}_{q}(\bm{x},p|\bm{x}_{0})=0$, one can employ the
eigenfunctions of $\mathcal{M}_{p}$ for a spectral decomposition of the
restriction of $\tilde{g}_{q}(\bm{x},p|\bm{x}_{0})$ to $\partial\Omega_{R}$,
which can then be extended to the whole domain $\Omega$. Finally, one takes
the inverse Laplace transform of $\tilde{G}_{q}(\bm{x},p|\bm{x}_{0})$ with
respect to $q$ to get the spectral expansion of the full propagator in the
Laplace domain:

$\displaystyle\tilde{P}(\bm{x},\ell,p|\bm{x}_{0})$
$\displaystyle=\tilde{G}_{\infty}(\bm{x},p|\bm{x}_{0})\delta(\ell)$ (20)
$\displaystyle+\frac{1}{D}\sum\limits_{k=0}^{\infty}[V_{k}^{(p)}(\bm{x}_{0})]^{*}\,V_{k}^{(p)}(\bm{x})\,e^{-\mu_{k}^{(p)}\ell},$

where asterisk denotes complex conjugate, and

$V_{k}^{(p)}(\bm{x}_{0})=\int\limits_{\partial\Omega_{R}}d\bm{x}\,\tilde{j}_{\infty}(\bm{x},p|\bm{x}_{0})\,v_{k}^{(p)}(\bm{x})$
(21)

is the extension of the eigenfunction $v_{k}^{(p)}(\bm{x})$ (defined on
$\partial\Omega_{R}$) to the whole domain $\Omega$. While the structure of the
spectral expansion (20) is identical to that derived in Grebenkov20 , its
“ingredients” $\tilde{G}_{\infty}(\bm{x},p|\bm{x}_{0})$,
$V_{k}^{(p)}(\bm{x}_{0})$ and $\mu_{k}^{(p)}$ are adapted to account for the
presence of reflecting and escape regions.

The inverse Laplace transform of Eq. (20) with respect to $p$ formally yields

$\displaystyle P(\bm{x},\ell,t|\bm{x}_{0})$
$\displaystyle=G_{\infty}(\bm{x},t|\bm{x}_{0})\delta(\ell)$ (22)
$\displaystyle+\mathcal{L}_{p\to
t}^{-1}\biggl{\\{}\frac{1}{D}\sum\limits_{k=0}^{\infty}[V_{k}^{(p)}(\bm{x}_{0})]^{*}\,V_{k}^{(p)}(\bm{x})\,e^{-\mu_{k}^{(p)}\ell}\biggr{\\}}.$

The first term in Eq. (22) represents the contribution of random trajectories
that moved from $\bm{x}_{0}$ to $\bm{x}$ in time $t$ without hitting neither
the target region $\partial\Omega_{R}$, nor the escape region
$\partial\Omega_{D}$ (in turn, they could encounter the reflecting part
$\partial\Omega_{N}$). Their “fraction” is precisely given by
$G_{\infty}(\bm{x},t|\bm{x}_{0})$. As the boundary local time $\ell_{t}^{R}$
remained zero for these trajectories, one gets the singular factor
$\delta(\ell)$. In turn, the second term accounts for all trajectories that
have encountered the target region $\partial\Omega_{R}$, but still avoided the
escape through $\partial\Omega_{D}$. The spectral representation (22) is an
alternative way for computing the full propagator, which is complementary to
Eq. (14). While both expressions involve an inverse Laplace transform, the
spectral characteristics of the Dirichlet-to-Neumann operator are in general
easier to access than the conventional propagator $G_{q}(\bm{x},t|\bm{x}_{0})$
in time domain. For instance, we will employ Eq. (22) in Sec. III to deduce
various properties of diffusion-controlled reactions in a spherical domain. We
emphasize that two representations provide complementary insights onto the
full propagator.

One can easily check that the functions $V_{k}^{(p)}(\bm{x})$ defined by Eq.
(21) satisfy:

$\displaystyle(p-D\Delta)V_{k}^{(p)}(\bm{x})$
$\displaystyle=0\quad(\bm{x}\in\Omega),$ (23a) $\displaystyle
V_{k}^{(p)}(\bm{x})$
$\displaystyle=v_{k}^{(p)}(\bm{x})\quad(\bm{x}\in\partial\Omega_{R}),$ (23b)
$\displaystyle V_{k}^{(p)}(\bm{x})$
$\displaystyle=0\quad(\bm{x}\in\partial\Omega_{D}),$ (23c)
$\displaystyle\partial_{n}V_{k}^{(p)}(\bm{x})$
$\displaystyle=0\quad(\bm{x}\in\partial\Omega_{N}).$ (23d)

Moreover, since $v_{k}^{(p)}$ is an eigenfunction of the Dirichlet-to-Neumann
operator, one has

$\partial_{n}V_{k}^{(p)}(\bm{x})=\mu_{k}^{(p)}V_{k}^{(p)}(\bm{x})=\mu_{k}^{(p)}v_{k}^{(p)}(\bm{x})\quad(\bm{x}\in\partial\Omega_{R}).$
(24)

Once the eigenfunction $v_{k}^{(p)}(\bm{x})$ is found, one can determine its
extension $V_{k}^{(p)}(\bm{x})$ either via Eq. (21), or by solving the above
problem (23). Alternatively, without knowing $v_{k}^{(p)}(\bm{x})$, one can
look directly at the eigenvalue problem (23a, 23c, 23d, 24), in which the
spectral parameter (here, $\mu_{k}^{(p)}$) stands in the boundary condition.
This is known as the Steklov problem (see Girouard17  and references therein),
while $\mu_{k}^{(p)}$ and $V_{k}^{(p)}(\bm{x})$ are the eigenvalues and
eigenfunctions of this problem. Despite this equivalence, we keep referring to
the Dirichlet-to-Neumann operator $\mathcal{M}_{p}$ and its spectral
properties.

The full propagator determines the corresponding probability flux density on
the boundary $\partial\Omega$:

$j(\bm{x},\ell,t|\bm{x}_{0})=-D\partial_{n}P(\bm{x},\ell,p|\bm{x}_{0})\quad(\bm{x}\in\partial\Omega).$
(25)

The spectral expansion (20) gives access to this quantity in the Laplace
domain

$\displaystyle\tilde{j}(\bm{x},\ell,p|\bm{x}_{0})$
$\displaystyle=\tilde{j}_{\infty}(\bm{x},p|\bm{x}_{0})\delta(\ell)$
$\displaystyle-\sum\limits_{k=0}^{\infty}[V_{k}^{(p)}(\bm{x}_{0})]^{*}\,(\partial_{n}V_{k}^{(p)})(\bm{x})\,e^{-\mu_{k}^{(p)}\ell}.$
(26)

The probability flux density $j(\bm{x},\ell,t|\bm{x}_{0})$ is the main object
of our study. In particular, we aim at providing its probabilistic
interpretation and deducing various related characteristics of diffusion-
controlled reactions in the presence of escape events. We will show that the
interpretation of $j(\bm{x},\ell,t|\bm{x}_{0})$ is more subtle than that of
$j_{q}(\bm{x},t|\bm{x}_{0})$ mentioned in Sec. II.1. For instance, it is easy
to prove that

$\int\limits_{0}^{\infty}d\ell\,j(\bm{x},\ell,t|\bm{x}_{0})=0\quad(\bm{x}\in\partial\Omega_{R}),$
(27)

i.e., the restriction of $j(\bm{x},\ell,t|\bm{x}_{0})$ to $\partial\Omega_{R}$
is not necessarily positive.

To show this relation, let us first integrate Eq. (20) over $\ell$ from $0$ to
infinity to get the following identity

$\displaystyle\tilde{G}_{0}(\bm{x},p|\bm{x}_{0})$
$\displaystyle=\tilde{G}_{\infty}(\bm{x},p|\bm{x}_{0})$ (28)
$\displaystyle+\frac{1}{D}\sum\limits_{k=0}^{\infty}[V_{k}^{(p)}(\bm{x}_{0})]^{*}\,V_{k}^{(p)}(\bm{x})\frac{1}{\mu_{k}^{(p)}}\,,$

so that

$\displaystyle\tilde{j}_{0}(\bm{x},p|\bm{x}_{0})$
$\displaystyle=\tilde{j}_{\infty}(\bm{x},p|\bm{x}_{0})$ (29)
$\displaystyle-\sum\limits_{k=0}^{\infty}[V_{k}^{(p)}(\bm{x}_{0})]^{*}\,(\partial_{n}V_{k}^{(p)})(\bm{x})\frac{1}{\mu_{k}^{(p)}}\,.$

For any $\bm{x}\in\partial\Omega_{R}$, the left-hand side of Eq. (29) is zero,
implying

$\tilde{j}_{\infty}(\bm{x},p|\bm{x}_{0})=\sum\limits_{k=0}^{\infty}[V_{k}^{(p)}(\bm{x}_{0})]^{*}\,v_{k}^{(p)}(\bm{x})\quad(\bm{x}\in\partial\Omega_{R}),$
(30)

where we applied Eq. (24). Using the identity (30), one can easily check that
the integral of Eq. (26) over $\ell$ from $0$ to infinity is strictly zero for
any $\bm{x}\in\partial\Omega_{R}$, that reads in time domain as Eq. (27).

In order to clarify the probabilistic meaning of
$j(\bm{x},\ell,t|\bm{x}_{0})$, we first look at the distribution of the
boundary local time $\ell_{t}^{R}$, then discuss the restriction of
$j(\bm{x},\ell,t|\bm{x}_{0})$ to the escape region $\partial\Omega_{D}$, and
finally describe its restriction to the target region $\partial\Omega_{R}$.

###  II.3 Probability density of the boundary local time

It is convenient to start by inspecting the distribution of the boundary local
time $\ell_{t}^{R}$. By definition, the integral of the full propagator
$P(\bm{x},\ell,t|\bm{x}_{0})$ over $\bm{x}\in\Omega$ determines the (marginal)
probability density of the boundary local time

$\rho(\ell,t|\bm{x}_{0})=\int\limits_{\Omega}d\bm{x}\,P(\bm{x},\ell,t|\bm{x}_{0}).$
(31)

Since the particle may escape the domain, this probability density is not
normalized to $1$:

$\int\limits_{0}^{\infty}d\ell\,\rho(\ell,t|\bm{x}_{0})=S_{0}(t|\bm{x}_{0})=\mathbb{P}_{\bm{x}_{0}}\\{T_{0}>t\\},$
(32)

i.e., the probability of no escape up to time $t$ (the subscript $0$ of
$T_{0}$ highlights that the target region $\partial\Omega_{R}$ is treated here
as reflecting). When there is no escape region, the survival probability
$S_{0}(t|\bm{x}_{0})$ is equal to $1$, ensuring the normalization in the
conventional case Grebenkov19a .

Integrating the fundamental relation (13) over $\bm{x}\in\Omega$, one gets

$S_{q}(t|\bm{x}_{0})=\int\limits_{0}^{\infty}d\ell\,e^{-q\ell}\,\rho(\ell,t|\bm{x}_{0})=\mathbb{E}_{\bm{x}_{0}}\\{e^{-q\ell_{t}^{R}}\\},$
(33)

i.e., the survival probability is the generating function of $\ell_{t}^{R}$.
Moreover, the inverse Laplace transform of Eq. (9) with respect to $q$ yields
another continuity equation

$\partial_{t}\rho(\ell,t|\bm{x}_{0})=-\bigl{[}J_{R}(\ell,t|\bm{x}_{0})+J_{D}(\ell,t|\bm{x}_{0})\bigr{]},$
(34)

where

$\displaystyle J_{R}(\ell,t|\bm{x}_{0})$
$\displaystyle=\int\limits_{\partial\Omega_{R}}d\bm{x}\,j(\bm{x},\ell,t|\bm{x}_{0}),$
(35a) $\displaystyle J_{D}(\ell,t|\bm{x}_{0})$
$\displaystyle=\int\limits_{\partial\Omega_{D}}d\bm{x}\,j(\bm{x},\ell,t|\bm{x}_{0}),$
(35b)

and we used that

$j_{q}(\bm{x},t|\bm{x}_{0})=\int\limits_{0}^{\infty}d\ell\,e^{-q\ell}\,j(\bm{x},\ell,t|\bm{x}_{0})$
(36)

due to Eq. (13).

In the Laplace domain, the integral of Eq. (20) yields

$\displaystyle\tilde{\rho}(\ell,p|\bm{x}_{0})$
$\displaystyle=\tilde{S}_{\infty}(p|\bm{x}_{0})\delta(\ell)$ (37)
$\displaystyle+\frac{1}{D}\sum\limits_{k=0}^{\infty}[V_{k}^{(p)}(\bm{x}_{0})]^{*}\,e^{-\mu_{k}^{(p)}\ell}\int\limits_{\Omega}d\bm{x}\,V_{k}^{(p)}(\bm{x}).$

###  II.4 Escape events

In a direct analogy with $j_{q}(\bm{x},t|\bm{x}_{0})$, the restriction of the
probability flux density $j(\bm{x},\ell,t|\bm{x}_{0})$ to the escape region
$\partial\Omega_{D}$ determines the joint probability density of the position
$\bm{X}_{T_{0}}$, the boundary local time $\ell_{T_{0}}^{R}$, and the escape
time $T_{0}$ defined by Eq. (5). If the exact location of the escape does not
matter, one can integrate over $\bm{x}\in\partial\Omega_{D}$ to get the joint
probability density $J_{D}(\ell,t|\bm{x}_{0})$ of $\ell_{T_{0}}^{R}$ and
$T_{0}$, see Eq. (35b). In the Laplace domain, this quantity reads

$\displaystyle\tilde{J}_{D}(\ell,p|\bm{x}_{0})$
$\displaystyle=\tilde{J}_{\infty}^{D}(p|\bm{x}_{0})\delta(\ell)$ (38)
$\displaystyle+\sum\limits_{k=0}^{\infty}[V_{k}^{(p)}(\bm{x}_{0})]^{*}\,e^{-\mu_{k}^{(p)}\ell}C_{k}^{(p)},$

where

$C_{k}^{(p)}=-\int\limits_{\partial\Omega_{D}}d\bm{x}\,\partial_{n}V_{k}^{(p)}(\bm{x})$
(39)

and $\tilde{J}_{\infty}^{D}(p|\bm{x}_{0})$ is given by Eq. (7a) at $q=\infty$.
While the escape time $T_{0}$ has been studied in the past, the joint
distribution of $\ell_{T_{0}}^{R}$ and $T_{0}$ has not been reported earlier.

Integrating Eq. (23a), using the Green’s formula and expressions (23d, 24),
one gets another representation:

$C_{k}^{(p)}=\mu_{k}^{(p)}\int\limits_{\partial\Omega_{R}}d\bm{x}\,v_{k}^{(p)}(\bm{x})-\frac{p}{D}\int\limits_{\Omega}d\bm{x}\,V_{k}^{(p)}(\bm{x}).$
(40)

One can apply this representation to check the correct normalization of
$J_{D}(\ell,t|\bm{x}_{0})$:

$\displaystyle\int\limits_{0}^{\infty}d\ell\int\limits_{0}^{\infty}dt\,J_{D}(\ell,t|\bm{x}_{0})=\int\limits_{0}^{\infty}d\ell\,\tilde{J}_{D}(\ell,p|\bm{x}_{0})$
$\displaystyle=\tilde{J}_{\infty}^{D}(0|\bm{x}_{0})+\sum\limits_{k=0}^{\infty}[V_{k}^{(0)}(\bm{x}_{0})]^{*}\frac{C_{k}^{(0)}}{\mu_{k}^{(0)}}$
$\displaystyle=\tilde{J}_{\infty}^{D}(0|\bm{x}_{0})+\tilde{J}_{\infty}^{R}(0|\bm{x}_{0})=1,$

where we used Eq. (30). The last equality follows from that the fact that
$\tilde{J}_{\infty}^{D}(0|\bm{x}_{0})$ and
$\tilde{J}_{\infty}^{R}(0|\bm{x}_{0})$ are the splitting probabilities (e.g.,
$\tilde{J}_{\infty}^{D}(0|\bm{x}_{0})$ is the probability of hitting
$\partial\Omega_{D}$ before hitting $\partial\Omega_{R}$) so that their sum is
equal to $1$.

The integral of $J_{D}(\ell,t|\bm{x}_{0})$ over $\ell$ determines the
(marginal) probability density function of the escape time $T_{0}$

$J_{0}^{D}(t|\bm{x}_{0})=\int\limits_{0}^{\infty}d\ell\,J_{D}(\ell,t|\bm{x}_{0}).$
(41)

Integrating the continuity relation (34) over $\ell$ from $0$ to infinity and
using the identity (27), one sees that

$\partial_{t}S_{0}(t|\bm{x}_{0})=-J_{0}^{D}(t|\bm{x}_{0}),$ (42)

where we also used the normalization (32). We thus re-derived that the
probability density of the first-passage time to $\partial\Omega_{D}$ is
obtained as the negative time derivative of the survival probability
$S_{0}(t|\bm{x}_{0})$. Expectedly, the derivation of this classical quantity
does not require the encounter-based approach.

In turn, the integral of $J_{D}(\ell,t|\bm{x}_{0})$ over $t$ yields the
(marginal) probability density of the boundary local time $\ell_{T_{0}}^{R}$
acquired before the escape:

$\displaystyle\rho_{D}(\ell|\bm{x}_{0})$
$\displaystyle=\int\limits_{0}^{\infty}dt\,J_{D}(\ell,t|\bm{x}_{0})=\tilde{J}_{D}(\ell,0|\bm{x}_{0})$
(43)
$\displaystyle=\tilde{J}_{\infty}^{D}(0|\bm{x}_{0})\delta(\ell)+\sum\limits_{k=0}^{\infty}[V_{k}^{(0)}(\bm{x}_{0})]^{*}\,e^{-\mu_{k}^{(0)}\ell}\,C_{k}^{(0)}.$

Here, the first term accounts for the trajectories that never hit the target
region $\partial\Omega_{R}$ and moved directly to the escape region (with the
splitting probability $\tilde{J}_{\infty}^{D}(0|\bm{x}_{0})$) so that the
acquired boundary local time is zero. In turn, the second term includes the
contribution of remaining trajectories.

The expression (38) allows one to determine joint positive integer-order
moments of $\ell_{T_{0}}^{R}$ and $T_{0}$ as

$\mathbb{E}_{\bm{x}_{0}}\\{[\ell_{T_{0}}^{R}]^{m}\,T_{0}^{n}\\}=(-1)^{n}m!\lim\limits_{p\to
0}\frac{\partial^{n}}{\partial
p^{n}}\sum\limits_{k=0}^{\infty}\frac{[V_{k}^{(p)}(\bm{x}_{0})]^{*}C_{k}^{(p)}}{[\mu_{k}^{(p)}]^{m+1}}$
(44)

for any $m=1,2,\ldots$ and any $n=0,1,\ldots$ (for $m=0$, the contribution
from the first (singular) term in Eq. (38) has to be included). In this way,
one can evaluate not only the moments of $\ell_{T_{0}}^{R}$ and $T_{0}$, but
also correlations between them. In particular, the mean of $\ell_{T_{0}}^{R}$
reads

$\displaystyle\mathbb{E}_{\bm{x}_{0}}\\{\ell_{T_{0}}^{R}\\}$
$\displaystyle=\sum\limits_{k=0}^{\infty}[V_{k}^{(0)}(\bm{x}_{0})]^{*}\frac{C_{k}^{(0)}}{[\mu_{k}^{(0)}]^{2}}$
$\displaystyle=\sum\limits_{k=0}^{\infty}[V_{k}^{(0)}(\bm{x}_{0})]^{*}\frac{1}{\mu_{k}^{(0)}}\int\limits_{\partial\Omega_{R}}d\bm{x}\,v_{k}^{(0)}(\bm{x})$
$\displaystyle=\int\limits_{\partial\Omega_{R}}d\bm{x}\,D\tilde{G}_{0}(\bm{x},0|\bm{x}_{0}),$
(45)

where we used the identity (28). This result is rather intuitive. In fact, the
Green’s function $D\tilde{G}_{0}(\bm{x},0|\bm{x}_{0})$ is known to be the mean
residence time of Brownian motion in a point $\bm{x}$ before the escape
Morters . The integral of this quantity over a thin layer of width $a$ near
the target region $\partial\Omega_{R}$ yields the mean residence time in this
layer, while the rescaling by $D/a$ results in the mean boundary local time in
the limit $a\to 0$, see Eq. (11).

To evaluate higher-order moments, we first integrate Eq. (20) $m+1$ times to
get

$\displaystyle\sum\limits_{k=0}^{\infty}\frac{[V_{k}^{(p)}(\bm{x}_{0})]^{*}\,V_{k}^{(p)}(\bm{x})}{[\mu_{k}^{(p)}]^{m+1}}$
$\displaystyle\quad=\int\limits_{0}^{\infty}d\ell_{1}\int\limits_{\ell_{1}}^{\infty}d\ell_{2}\,\ldots\int\limits_{\ell_{m+1}}^{\infty}d\ell_{m+1}D\tilde{P}(\bm{x},\ell_{m+1},p|\bm{x}_{0})$
$\displaystyle\quad=\frac{1}{m!}\int\limits_{0}^{\infty}d\ell\,\ell^{m}\,D\tilde{P}(\bm{x},\ell,p|\bm{x}_{0})$
(46)

for any $m=1,2,\ldots$. In turn, differentiating $m$ times the fundamental
relation (13) with respect to $q$, we derive the following identity:

$\sum\limits_{k=0}^{\infty}\frac{[V_{k}^{(p)}(\bm{x}_{0})]^{*}\,V_{k}^{(p)}(\bm{x})}{[\mu_{k}^{(p)}]^{m+1}}=\frac{(-1)^{m}}{m!}\lim\limits_{q\to
0}\partial_{q}^{m}D\tilde{G}_{q}(\bm{x},p|\bm{x}_{0})$ (47)

for $m=1,2,\ldots$ (for $m=0$, there is an additional term, see Eq. (28)).
This identity provides a peculiar interpretation of the derivatives of the
Green’s function $\tilde{G}_{q}(\bm{x},p|\bm{x}_{0})$ with respect to $q$.
Using Eqs. (39, 44), one gets

$\displaystyle\mathbb{E}_{\bm{x}_{0}}\\{[\ell_{T_{0}}^{R}]^{m}\,T_{0}^{n}\\}$
$\displaystyle=\frac{(-1)^{m+n+1}}{m!}$ (48)
$\displaystyle\times\lim\limits_{q\to 0\atop p\to
0}\partial^{m}_{q}\partial^{n}_{p}\int\limits_{\partial\Omega_{D}}d\bm{x}\,D\partial_{n}\tilde{G}_{q}(\bm{x},p|\bm{x}_{0}).$

The definitions (6, 7a) further simplify this relation as

$\mathbb{E}_{\bm{x}_{0}}\\{[\ell_{T_{0}}^{R}]^{m}\,T_{0}^{n}\\}=(-1)^{m+n}\lim\limits_{q\to
0\atop p\to
0}\partial^{m}_{q}\partial^{n}_{p}\tilde{J}_{q}^{D}(p|\bm{x}_{0}).$ (49)

Note that this relation is also applicable for $m=0$, in which case it reduces
to the standard expression for the moments of the first-passage time $T_{0}$:

$\mathbb{E}_{\bm{x}_{0}}\\{T_{0}^{n}\\}=(-1)^{n}\lim\limits_{p\to
0}\partial^{n}_{p}\tilde{J}_{0}^{D}(p|\bm{x}_{0}).$ (50)

Like $\tilde{J}_{0}^{D}(p|\bm{x}_{0})$ was known to be the generating function
of the first-passage time $T_{0}$, $\tilde{J}_{q}^{D}(p|\bm{x}_{0})$ turns out
to be the generating function of both $T_{0}$ and $\ell_{T_{0}}^{R}$. This is
not surprising given that $p$ and $q$ appear as the conjugate variables in the
Laplace transforms with respect to $t$ and $\ell$. The elegant relation (49)
highlights the similarity between the physical time $t$, which can be seen as
a proxy for the number of “elementary jumps” of the particle in the bulk, and
the boundary local time $\ell$, which is a proxy of the number of its “jumps”
on the target region $\partial\Omega_{R}$. This similarity was already
discussed in Grebenkov20  and particularly in Grebenkov20b , in which a
surface-hopping propagator, based on the boundary local time, was introduced.

###  II.5 Reaction events

Now we can inspect the restriction of the probability flux density
$j(\bm{x},\ell,t|\bm{x}_{0})$ to the target region $\partial\Omega_{R}$. In
the Laplace domain, the spectral expansion (20) implies

$\displaystyle\tilde{j}(\bm{x},\ell,p|\bm{x}_{0})$
$\displaystyle=\tilde{j}_{\infty}(\bm{x},p|\bm{x}_{0})\delta(\ell)$ (51)
$\displaystyle-\sum\limits_{k=0}^{\infty}[V_{k}^{(p)}(\bm{x}_{0})]^{*}\,\mu_{k}^{(p)}v_{k}^{(p)}(\bm{x})e^{-\mu_{k}^{(p)}\ell}$

for any $\bm{x}\in\partial\Omega_{R}$. While this relation allows one to
compute this quantity, its probabilistic interpretation remains tricky, in
particular, due to its negative values according to Eq. (27).

To clarify this point, we first rewrite Eq. (4a) as
$DG_{q}(\bm{x},t|\bm{x}_{0})=j_{q}(\bm{x},t|\bm{x}_{0})/q$ on
$\partial\Omega_{R}$ and take its inverse Laplace transform with respect to
$q$ to get

$DP(\bm{x},\ell,t|\bm{x}_{0})=\int\limits_{0}^{\ell}d\ell^{\prime}\,j(\bm{x},\ell^{\prime},t|\bm{x}_{0})\quad(\bm{x}\in\partial\Omega_{R}).$
(52)

The probability flux density $j(\bm{x},\ell^{\prime},t|\bm{x}_{0})$ on
$\partial\Omega_{R}$ can thus also be seen as the derivative of the full
propagator on $\partial\Omega_{R}$ with respect to $\ell$. This is consistent
with the boundary value problem for the full propagator discussed in
Grebenkov20 . In the limit $\ell\to\infty$, Eq. (52) reduces to the relation
(27) because $P(\bm{x},\ell,t|\bm{x}_{0})\to 0$, i.e., the probability of
getting infinitely large values of the boundary local time $\ell_{t}^{R}$ at a
finite time $t$ is zero. As the left-hand side of Eq. (52) is a probability
density, the integral in the right-hand side is nonnegative, despite eventual
negative values of $j(\bm{x},\ell^{\prime},t|\bm{x}_{0})$. What does it
represent?

In order to get its probabilistic interpretation, we introduce the first-
crossing time $\mathcal{T}_{\ell}$ of a threshold $\ell$ by the boundary local
time $\ell_{t}^{R}$:

$\mathcal{T}_{\ell}=\inf\\{t>0~{}:~{}\ell_{t}^{R}>\ell\\}.$ (53)

If the particle has escaped the domain before crossing the threshold, the
first-crossing time is set to infinity. We have then

$\mathbb{P}\\{\mathcal{T}_{\ell}>t\\}=\mathbb{P}\\{\ell_{t}^{R}<\ell,~{}T_{0}>t\\}+\mathbb{P}\\{\ell_{T_{0}}^{R}<\ell,~{}T_{0}<t\\}.$
(54)

The first term describes no crossing of the threshold $\ell$ when the escape
occurs after time $t$. In turn, the second term describes the escape event
before $t$, for which the acquired boundary local time $\ell_{T_{0}}^{R}$
remains below the threshold. The probability density of the first-crossing
time reads then as

$\displaystyle U(\ell,t|\bm{x}_{0})$
$\displaystyle=-\partial_{t}\mathbb{P}\\{\mathcal{T}_{\ell}>t\\}$
$\displaystyle=-\partial_{t}\int\limits_{0}^{\ell}d\ell^{\prime}\,\rho(\ell^{\prime},t|\bm{x}_{0})-\int\limits_{0}^{\ell}d\ell^{\prime}\,J_{D}(\ell^{\prime},t|\bm{x}_{0}),$

where the probability density $\rho(\ell^{\prime},t|\bm{x}_{0})$ was used to
evaluate the first term in Eq. (54), and the joint probability density
$J_{D}(\ell,t|\bm{x}_{0})$ of $\ell_{T_{0}}^{R}$ and $T_{0}$ for the second
term. Note that the definition of $\rho(\ell^{\prime},t|\bm{x}_{0})$
automatically accounts for no escape up to time $t$. Using the continuity
equation (34), we finally get

$U(\ell,t|\bm{x}_{0})=\int\limits_{0}^{\ell}d\ell^{\prime}\,J_{R}(\ell^{\prime},t|\bm{x}_{0}).$
(55)

Despite the fact that $J_{R}(\ell^{\prime},t|\bm{x}_{0})$ cannot be
interpreted as a joint probability density (in analogy with
$J_{D}(\ell,t|\bm{x}_{0})$), its integral over $\ell^{\prime}$ yields the
probability density of the first-crossing time $\mathcal{T}_{\ell}$.

According to Eq. (52), one gets another representation

$U(\ell,t|\bm{x}_{0})=\int\limits_{\partial\Omega_{R}}d\bm{x}\,DP(\bm{x},\ell,t|\bm{x}_{0}),$
(56)

which was earlier derived in Grebenkov20  for the particular case
$\partial\Omega_{R}=\partial\Omega$ (i.e., without the escape region). This
relation has an intuitive interpretation. Let us again consider a thin layer
of width $a$ near $\partial\Omega_{R}$. By definition, the integral of the
full propagator $P(\bm{x},\ell,t|\bm{x}_{0})$ over
$\bm{x}\in\partial\Omega_{R}$, multiplied by $a$ and $d\ell$, is the
probability of finding the particle in that layer at time $t$ with the
boundary local time $\ell_{t}^{R}$ belonging to $(\ell,\ell+d\ell)$. As
$\ell_{t}^{R}$ is an nondecreasing process that increments only when
$\bm{X}_{t}\in\partial\Omega_{R}$, the value $\ell$ is thus achieved for the
first time at $t$, i.e.,

$\displaystyle
a\,d\ell\int\limits_{\partial\Omega_{R}}d\bm{x}\,P(\bm{x},\ell,t|\bm{x}_{0})$
$\displaystyle=\mathbb{P}_{\bm{x}_{0}}\\{\mathcal{T}_{\ell}\in(t,t+dt)\\}$
$\displaystyle=U(t,\ell|\bm{x}_{0})dt.$ (57)

Since the increments $d\ell$ and $dt$ are related as $d\ell=Ddt/a$ according
to Eq. (12), one gets Eq. (56).

In the Laplace domain, Eq. (51) implies

$\displaystyle\tilde{U}(\ell,p|\bm{x}_{0})$
$\displaystyle=\sum\limits_{k=0}^{\infty}[V_{k}^{(p)}(\bm{x}_{0})]^{*}\,e^{-\mu_{k}^{(p)}\ell}\int\limits_{\partial\Omega_{R}}d\bm{x}\,v_{k}^{(p)}(\bm{x}).$

Setting $\ell=0$ and employing Eq. (30), one gets

$\displaystyle\tilde{U}(0,p|\bm{x}_{0})$
$\displaystyle=\sum\limits_{k=0}^{\infty}[V_{k}^{(p)}(\bm{x}_{0})]^{*}\,\int\limits_{\partial\Omega_{R}}d\bm{x}\,v_{k}^{(p)}(\bm{x})$
$\displaystyle=\int\limits_{\partial\Omega_{R}}d\bm{x}\,\tilde{j}_{\infty}(\bm{x},p|\bm{x}_{0})=\tilde{J}_{\infty}^{R}(p|\bm{x}_{0}),$

i.e.,

$U(0,t|\bm{x}_{0})=J_{\infty}^{R}(t|\bm{x}_{0}).$ (58)

In other words, as the first crossing of the threshold $\ell=0$ corresponds to
the first arrival onto $\partial\Omega_{R}$, the first-crossing time
$\mathcal{T}_{0}$ is simply the first-passage time to $\partial\Omega_{R}$
(before escaping the domain).

As earlier stressed in Grebenkov20 , the probability density
$U(\ell,t|\bm{x}_{0})$ of the first-crossing time is tightly related to the
probability density $J_{q}^{R}(t|\bm{x}_{0})$ of the reaction time $\tau_{q}$
on $\partial\Omega_{R}$ discussed in Sec. II.1. In fact, using Eqs. (7b, 35a,
36), we first get

$J_{q}^{R}(t|\bm{x}_{0})=\int\limits_{0}^{\infty}d\ell\,e^{-q\ell}\,J_{R}(\ell,t|\bm{x}_{0}).$
(59)

In order to transform $J_{R}(\ell,t|\bm{x}_{0})$ into $U(\ell,t|\bm{x}_{0})$,
one needs to integrate by parts. However, $J_{R}(\ell,t|\bm{x}_{0})$ contains
a singular term $J_{\infty}^{R}(t|\bm{x}_{0})\delta(\ell)$ that has to be
treated separately. We have then

$\displaystyle J_{q}^{R}(t|\bm{x}_{0})$
$\displaystyle=J_{\infty}^{R}(t|\bm{x}_{0})+\int\limits_{0^{+}}^{\infty}d\ell\,e^{-q\ell}\,J_{R}(\ell,t|\bm{x}_{0})$
$\displaystyle=J_{\infty}^{R}(t|\bm{x}_{0})+\underbrace{\left.\biggl{(}\int\limits_{0^{+}}^{\infty}d\ell\,e^{-q\ell}\,U(\ell,t|\bm{x}_{0})\biggr{)}\right|_{\ell=0}^{\infty}}_{-U(0,t|\bm{x}_{0})}$
$\displaystyle+q\int\limits_{0^{+}}^{\infty}d\ell\,e^{-q\ell}\,U(\ell,t|\bm{x}_{0}),$

where we wrote the lower bound of the integral as $0^{+}$ to highlight that
the singular term was excluded. According to Eq. (58), two first terms cancel
each other, yielding

$J_{q}^{R}(t|\bm{x}_{0})=\int\limits_{0}^{\infty}d\ell\,q\,e^{-q\ell}\,U(\ell,t|\bm{x}_{0}).$
(60)

To interpret this relation, we introduce a random threshold $\hat{\ell}$
obeying the exponential probability law with the rate $q$,
$\mathbb{P}\\{\hat{\ell}>\ell\\}=e^{-q\ell}$, so that $qe^{-q\ell}$ is the
probability density of this law. As a consequence, the integral (60) is the
average over random realizations of the threshold $\hat{\ell}$ of the
probability density of the first-crossing time $\mathcal{T}_{\hat{\ell}}$.
Since the left-hand side is the probability density of the reaction time
$\tau_{q}$, we conclude that

$\tau_{q}=\mathcal{T}_{\hat{\ell}}=\inf\\{t>0~{}:~{}\ell_{t}^{R}>\hat{\ell}\\},$
(61)

i.e., the reaction occurs when the boundary local time $\ell_{t}^{R}$ exceeds
the random threshold $\hat{\ell}$ with the exponential law. We emphasize that
the exponential law follows directly from the postulated Robin boundary
condition (4a). This interpretation, earlier suggested in Grebenkov20 , allows
one to go beyond the Robin boundary condition and to implement various surface
reaction mechanisms by choosing an appropriate law for the random threshold
$\hat{\ell}$. In this way, our definition (61) of the reaction time $\tau_{q}$
remains valid in a much more general setting, while Eq. (60), which was
specific to the Robin boundary condition, is then generalized to

$J_{\psi}^{R}(t|\bm{x}_{0})=\int\limits_{0}^{\infty}d\ell\,\psi(\ell)\,U(\ell,t|\bm{x}_{0}),$
(62)

where $\psi(\ell)$ is the chosen probability density of the random threshold
$\hat{\ell}$. Different choices of $\psi(\ell)$ and its consequences on the
distribution of the reaction time were discussed in Grebenkov20 . This
extension of surface reaction mechanisms is directly applicable to our setting
with the escape region $\partial\Omega_{D}$. In this way, we made a further
step towards more realistic modeling of diffusion-controlled reactions by
incorporating the effect of escape events.

We complete this section by providing a deeper probabilistic interpretation of
the first-passage time $T_{q}$ to the escape region $\partial\Omega_{D}$ in
the presence of the reactive region $\partial\Omega_{R}$. As discussed in Sec.
II.1, this random variable is described by the probability density
$J_{q}^{D}(t|\bm{x}_{0})$, which is related to $J_{D}(\ell,t|\bm{x}_{0})$ due
to Eq. (36) as

$J_{q}^{D}(t|\bm{x}_{0})=\int\limits_{0}^{\infty}d\ell\,e^{-q\ell}\,J_{D}(\ell,t|\bm{x}_{0}).$
(63)

As previously, the factor $e^{-q\ell}$ incorporates the condition of no
reaction until the escape, which was automatically included into the left-hand
side via the Robin boundary condition on $\partial\Omega_{R}$. In turn, if the
reaction happens before the escape, the escape time $T_{q}$ is set to
infinity. In other words, one has

$T_{q}=\left\\{\begin{array}[]{ll}T_{0}&\textrm{if}~{}\ell_{T_{0}}^{R}<\hat{\ell},\\\
+\infty&\textrm{otherwise.}\end{array}\right.$ (64)

As earlier, the dependence on the reactivity parameter $q$ is incorporated via
the random threshold $\hat{\ell}$ obeying the exponential law with the rate
$q$. This representation highlights the effect of surface reactions onto the
fisrt-passage time $T_{q}$. In fact, if the target region $\partial\Omega_{R}$
was inert, the particle would reach the escape region $\partial\Omega_{D}$ at
a random time $T_{0}$. In turn, the reactivity of $\partial\Omega_{R}$ makes
the only change that the particle can react on $\partial\Omega_{R}$ and thus
never escape. If the boundary local time $\ell_{T_{0}}^{R}$ at $T_{0}$ has not
crossed the random threshold $\hat{\ell}$, the reaction on
$\partial\Omega_{R}$ has not happened, and $T_{q}=T_{0}$. In turn, if the
threshold $\hat{\ell}$ has been crossed, the reaction occurred before $T_{0}$
so that $T_{q}=\infty$. The above definition of the first-passage time to the
escape region $\partial\Omega_{D}$ within the encounter-based approach allows
one to study this quantity for other surface reaction mechanisms, beyond the
conventional one described by Robin boundary condition. For this purpose, one
can choose an appropriate law $\mathbb{P}\\{\hat{\ell}>\ell\\}$ for the random
threshold $\hat{\ell}$ and replace the factor $e^{-q\ell}$ in Eq. (63) by this
law. Further investigations of this setting present an interesting perspective
of the present work.

##  III Spherical target

In order to illustrate the general properties of diffusion-controlled
reactions with eventual escape, we consider diffusion between two concentric
spheres of radii $R$ and $L$:
$\Omega=\\{\bm{x}\in\mathbb{R}^{3}~{}:~{}R<|\bm{x}|<L\\}$. The inner sphere
represents the target region $\partial\Omega_{R}$, while the outer sphere is
the escape region $\partial\Omega_{D}$ (in this setting,
$\partial\Omega_{N}=\emptyset$). In spherical coordinates
$\bm{x}=(r,\theta,\phi)$, the modified Helmholtz equation can be solved via
separation of variables. In particular, the eigenbasis of the Dirichlet-to-
Neumann operator is well known Grebenkov20b ; Grebenkov20c ; Grebenkov22c . In
fact, the rotational invariance of the domain implies that the eigenfunctions
of $\mathcal{M}_{p}$ are the (normalized) spherical harmonics
$Y_{mn}(\theta,\phi)$

$v_{nm}=\frac{1}{R}Y_{mn}(\theta,\phi)\qquad\left({n=0,1,\ldots,\atop|m|\leq
n}\right),$ (65)

where the double index $nm$ is used instead of a single index $k$. In turn,
the eigenvalues are

$\mu_{nm}^{(p)}=-\bigl{(}\partial_{r}g_{n}^{(p)}(r)\bigr{)}_{r=R}\,,$ (66)

where

$g_{n}^{(p)}(r)=\frac{k_{n}(\alpha L)i_{n}(\alpha r)-i_{n}(\alpha
L)k_{n}(\alpha r)}{k_{n}(\alpha L)i_{n}(\alpha R)-i_{n}(\alpha L)k_{n}(\alpha
R)}$ (67)

are the radial functions, with $\alpha=\sqrt{p/D}$ and $i_{n}(z)$ and
$k_{n}(z)$ being the modified spherical Bessel functions of the first and
second kind. Note that $g_{n}^{(p)}(L)=0$ and $g_{n}^{(p)}(R)=1$. One sees
that the eigenvalues are $(2n+1)$ times degenerate (they do not depend on the
index $m$), while the eigenfunctions do not depend on $p$. Finally, it is easy
to check that

$V_{nm}^{(p)}(\bm{x})=g_{n}^{(p)}(r)\,v_{nm}(\theta,\phi).$ (68)

As most quantities of interest are obtained by integrating over the target
region $\partial\Omega_{R}$, the orthogonality of eigenfunctions $v_{nm}$ to
$v_{00}=1/\sqrt{4\pi R^{2}}$ ensures that all terms vanish except this ground
eigenmode. Since $i_{0}(z)=\sinh(z)/z$ and $k_{0}(z)=e^{-z}/z$, one finds

$g_{0}^{(p)}(r)=\frac{R\,\sinh(\alpha(L-r))}{r\,\sinh(\alpha(L-R))}$ (69)

and

$\mu_{0}^{(p)}=\mu_{00}^{(p)}=\frac{1}{R}+\alpha\,\mathrm{ctanh}(\alpha(L-R)).$
(70)

Note that

$\displaystyle\int\limits_{\Omega}d\bm{x}\,V_{00}^{(p)}(\bm{x})=\frac{4\pi}{\sqrt{4\pi
R^{2}}}\int\limits_{R}^{L}dr\,r^{2}\,g_{0}^{(p)}(r)$
$\displaystyle=\frac{\sqrt{4\pi}}{\sinh(\alpha(L-R))}\biggl{(}-\frac{L}{\alpha}$
$\displaystyle+\frac{\alpha
R\cosh(\alpha(L-R))+\sinh(\alpha(L-R))}{\alpha^{2}}\biggr{)}$

and

$C_{0}^{(p)}=\sqrt{4\pi}\,\frac{\alpha L}{\sinh(\alpha(L-R))}\,.$

One can also compute

$\displaystyle J_{q}^{D}(p|\bm{x}_{0})$
$\displaystyle=\frac{L}{r_{0}}\,\frac{(1+qR)\sinh(\alpha(r_{0}-R))+\alpha
R\cosh(\alpha(r_{0}-R))}{(1+qR)\sinh(\alpha(L-R))+\alpha
R\cosh(\alpha(L-R))}\,,$ (71a) $\displaystyle J_{q}^{R}(p|\bm{x}_{0})$
$\displaystyle=\frac{R}{r_{0}}\,\frac{qR\sinh(\alpha(L-r_{0}))}{(1+qR)\sinh(\alpha(L-R))+\alpha
R\cosh(\alpha(L-R))}\,,$ (71b)

where $r_{0}=|\bm{x}_{0}|$, from which

$\displaystyle J_{\infty}^{R}(p|\bm{x}_{0})$
$\displaystyle=\frac{R\,\sinh(\alpha(L-r_{0}))}{r_{0}\,\sinh(\alpha(L-R))}\,,$
$\displaystyle J_{\infty}^{D}(p|\bm{x}_{0})$
$\displaystyle=\frac{L\sinh(\alpha(r_{0}-R))}{r_{0}\sinh(\alpha(L-R))}\,,$
$\displaystyle J_{0}^{D}(p|\bm{x}_{0})$
$\displaystyle=\frac{L}{r_{0}}\,\frac{\sinh(\alpha(r_{0}-R))+\alpha
R\cosh(\alpha(r_{0}-R))}{\sinh(\alpha(L-R))+\alpha R\cosh(\alpha(L-R))}\,.$

In the following illustrations, we fix units of length and time by setting
$R=1$ and $D=1$.

###  III.1 Escape events

According to Eq. (38), we find

$\displaystyle\tilde{J}_{D}(\ell,p|\bm{x}_{0})$
$\displaystyle=\frac{L\sinh(\alpha(r_{0}-R))}{r_{0}\sinh(\alpha(L-R))}\delta(\ell)$
$\displaystyle+\frac{L\,\alpha\,\sinh(\alpha(L-r_{0}))}{r_{0}\,\sinh^{2}(\alpha(L-R))}\,e^{-\mu_{0}^{(p)}\ell}\,,$
(73)

from which an inverse Laplace transform with respect to $p$ yields the joint
probability density $J_{D}(\ell,t|\bm{x}_{0})$ of $\ell_{T_{0}}^{R}$ and
$T_{0}$. The inversion of the first term in Eq. (73) can be found explicitly
via the residue theorem. In turn, the presence of the $p$-dependent factor
$e^{-\mu_{0}^{(p)}\ell}$ in the second term makes the inversion challenging
(see Grebenkov20c  for some analytical tools). For this reason, we employed a
numerical Laplace transform inversion by the Talbot algorithm Talbot79 .

Figure 2 illustrates the behavior of $J_{D}(\ell,t|\bm{x}_{0})$. Two terms in
Eq. (73) give respectively singular and regular contributions to this density
with respect to $\ell$. The regular contribution as a function of $\ell$ and
$t$ is shown by a surface. In turn, the singular contribution corresponding to
$\ell=0$ is shown by a black line. For the panel (a), we set $L=2$ and
$r_{0}=1.5$, i.e., the particle starts in the middle between the target region
at $R=1$ and the escape region at $L=2$; in turn, in the panel (b), the escape
region is moved to $L=10$.

Figure 2:  Joint probability density $J_{D}(\ell,t|\bm{x}_{0})$ of the
boundary local time $\ell_{T_{0}}^{R}$ and the escape time $T_{0}$ for
diffusion between two concentric spheres of radii $R$ and $L$, with $R=1$,
$r_{0}=1.5$, $D=1$, $L=2$ (a) and $L=10$ (b). Surface shows the regular part
of this density, while black solid line presents the prefactor
$J_{\infty}^{D}(t|\bm{x}_{0})$ in front of the singular term $\delta(\ell)$.
Note that this curve should be located at $\ell=0$, which is not visible on
the logarithmic scale; it was thus artificially put at $\ell=10^{-2}$ for
illustration purposes. The ranges of times $t$ differ by factor $10$ between
two panels.

Since the particle needs to diffuse from its starting point $\bm{x}_{0}$ to
the escape region at $L$, its escape at very short times is highly unlikely.
In particular, the probability density $J_{\infty}^{D}(t|\bm{x}_{0})$ of the
first-passage time to $\partial\Omega_{D}$ determining the singular term, is
known to exhibit the short-time behavior of a Lévy-Smirnov form
$t^{-3/2}e^{-(L-r_{0})^{2}/(4Dt)}$ (see Godec16 ; Grebenkov18  and references
therein). Similar behavior can be derived for the regular contribution. In the
same vein, the probability of escape at very long times is also negligible
because the particle cannot avoid the escape region too long. These arguments
explain a distinct maximum of $J_{D}(\ell,t|\bm{x}_{0})$ with respect to $t$
at intermediate times. Finally, too large values of the boundary local time
$\ell_{T_{0}}^{R}$ for intermediate escape times would require for the
particle to stay too long near the target region $\partial\Omega_{R}$, which
is also highly improbable. One sees therefore that
$J_{\infty}^{D}(t|\bm{x}_{0})$ exhibits a single “boss” that rapidly goes down
as $t\to 0$, $t\to\infty$, and $\ell\to\infty$.

The panel (a) of Fig. 2 indicates that the maximum of the singular
contribution (the curve $J_{\infty}^{D}(t|\bm{x}_{0})$) is shifted to shorter
times with respect to the boss of $J_{D}(\ell,t|\bm{x}_{0})$. In fact, the
particle that never hit the target region $\partial\Omega_{R}$ (and thus has
$\ell_{T_{0}}^{R}=0$) needs to travel the distance $L-r_{0}$ from the starting
point to the escape region that determines the maximum of
$J_{\infty}^{D}(t|\bm{x}_{0})$. In turn, the particle that first encountered
$\partial\Omega_{R}$ and then escaped, has to travel the longer distance
$(r_{0}-R)+(L-R)$ that requires longer time. When the starting point gets
closer to $\partial\Omega_{D}$, the separation between two maxima is even
larger, while the location of $\bm{x}_{0}$ near $\partial\Omega_{R}$ would
reduce this separation. Similarly, if $L$ is much larger than $r_{0}$ (as on
the panel (b)), the difference between two distances is relatively small, and
two maxima are close.

Setting $p=0$ in Eq. (73), one gets the (marginal) probability density of
$\ell_{T_{0}}^{R}$:

$\displaystyle\rho_{D}(\ell|\bm{x}_{0})$
$\displaystyle=\frac{L}{r_{0}(L-R)}\biggl{(}(r_{0}-R)\delta(\ell)$
$\displaystyle+\frac{L-r_{0}}{L-R}\,e^{-\ell(1/R+1/(L-R))}\biggr{)}.$ (74)

This is a mixture of an exponential distribution and an atom at $\ell=0$.
Figure 3 illustrates the regular (exponential) part of this distribution for
the example from Fig. 2(a). As an additional validation step, we also realized
Monte Carlo simulations of random trajectories and thus independently
evaluated the statistics of $\ell_{T_{0}}^{R}$. For this purpose, a random
trajectory of a particle was simulated by adding independent Gaussian
increments, with mean zero and variance $\sqrt{2D\delta}$ in each direction,
where $\delta$ is a small time step. If the particle jumps inside the inner
sphere of radius $R$, it is normally reflected back. At each step when the
particle is within a layer of width $a$ near the inner sphere (i.e., when
$|\bm{X}_{t}|<R+a$), the boundary local time $\ell_{t}^{R}$ is incremented by
$D\delta/a$ according to Eq. (12). The simulation is stopped when the particle
crosses the outer sphere of radius $L$. Repeating such a simulation $M$ times,
one gets an empirical statistics of the escape time $T_{0}$ and of the
acquires boundary local time $\ell_{T_{0}}^{R}$.

Figure 3 confirms an excellent agreement between the exact solution (74) and
simulations. One can also appreciate that getting the behavior of
$\rho_{D}(\ell|\bm{x}_{0})$ at too small or too large $\ell$ is problematic
for Monte Carlo simulations. In fact, estimations at small $\ell$ require very
short time steps in the modeling of the random trajectory (and thus too long
simulations). In turn, estimations at large $\ell$ require too many simulated
trajectories, as the probability of getting large $\ell_{T_{0}}^{R}$ is
exponentially small.

Figure 3:  The regular part of the probability density
$\rho_{D}(\ell|\bm{x}_{0})$ of the acquired boundary local time
$\ell_{T_{0}}^{R}$ up to the escape time $T_{0}$ for diffusion between two
concentric spheres of radii $R$ and $L$, with $R=1$, $L=2$, $r_{0}=1.5$, and
$D=1$. Solid line shows the exact solution (74), while symbols present a
renormalized histogram from Monte Carlo simulations with $M=10^{4}$ particles,
the time step $\delta=10^{-4}$ and the layer width $a=5\sqrt{2D\delta}$.

The moments of $\ell_{T_{0}}^{R}$ are particularly simple:

$\mathbb{E}_{\bm{x}_{0}}\\{[\ell_{T_{0}}^{R}]^{m}\\}=m!\,\frac{R(L-r_{0})}{r_{0}(L-R)}\biggl{(}\frac{1}{R}+\frac{1}{L-R}\biggr{)}^{-m}.$
(75)

In particular,

$\displaystyle\mathbb{E}_{\bm{x}_{0}}\\{\ell_{T_{0}}^{R}\\}$
$\displaystyle=R^{2}(1/r_{0}-1/L)\,,$ (76a)
$\displaystyle\mathbb{V}_{\bm{x}_{0}}\\{\ell_{T_{0}}^{R}\\}$
$\displaystyle=R^{3}(1/r_{0}-1/L)(2-R/r_{0}-R/L),$ (76b)

where $\mathbb{V}_{\bm{x}_{0}}$ denotes the variance. The moments of $T_{0}$
are determined via Eq. (50); in particular, one has

$\mathbb{E}_{\bm{x}_{0}}\\{T_{0}\\}=\frac{(L-r_{0})\bigl{(}r_{0}L(r_{0}+L)-2R^{3}\bigr{)}}{6Dr_{0}L}$
(77)

and

$\displaystyle\mathbb{V}_{\bm{x}_{0}}\\{T_{0}\\}=\frac{L-r_{0}}{90D^{2}L^{2}r_{0}^{2}}\biggl{(}r_{0}^{2}L^{5}+L^{4}r_{0}^{3}+r_{0}^{4}L^{3}+r_{0}^{5}L^{2}$
$\displaystyle\quad-20L^{2}r_{0}^{2}R^{3}+36R^{5}r_{0}L-10R^{6}L-10R^{6}r_{0}\biggr{)}.$
(78)

In turn, the joint moments can be found from

$\displaystyle\mathbb{E}_{\bm{x}_{0}}\\{[\ell_{T_{0}}^{R}]^{m}\,T_{0}^{n}\\}=(-1)^{n}m!\frac{L}{r_{0}}$
$\displaystyle\quad\times\lim\limits_{p\to 0}\frac{\partial^{n}}{\partial
p^{n}}\biggl{(}\frac{\alpha\,\sinh(\alpha(L-r_{0}))}{\sinh^{2}(\alpha(L-R))}[\mu_{0}^{(p)}]^{-m-1}\biggr{)}.$
(79)

For instance, we get

$\displaystyle\mathbb{E}_{\bm{x}_{0}}\\{\ell_{T_{0}}^{R}\,T_{0}\\}$
$\displaystyle=\frac{R^{2}(L-r_{0})}{6DL^{2}r_{0}}\biggl{(}2(L+2R)(L-R)^{2}$
$\displaystyle-L(L-r_{0})^{2}\biggr{)}\,.$ (80)

Expectedly, all these quantities vanish as $r_{0}\to L$ because the particle
started on $\partial\Omega_{D}$ escapes immediately, yielding
$T_{0}=\ell_{T_{0}}^{R}=0$ in a deterministic way. The above expressions allow
one to compute the Peason’s correlation coefficient between $T_{0}$ and
$\ell_{T_{0}}^{R}$:

$C=\frac{\mathbb{E}_{\bm{x}_{0}}\\{\ell_{T_{0}}^{R}\,T_{0}\\}-\mathbb{E}_{\bm{x}_{0}}\\{\ell_{T_{0}}^{R}\\}\,\mathbb{E}_{\bm{x}_{0}}\\{T_{0}\\}}{\sqrt{\mathbb{V}_{\bm{x}_{0}}\\{T_{0}\\}}\,\sqrt{\mathbb{V}_{\bm{x}_{0}}\\{\ell_{T_{0}}^{R}\\}}}\,.$
(81)

Figure 4 shows the correlation coefficient $C$ as a function of $r_{0}/L$ for
two values of $L$. In both cases, the correlation is positive. Indeed, if the
particle escapes faster (at smaller $T_{0}$), it would generally have lower
chances to encounter the target region frequently so that the boundary local
time $\ell_{T_{0}}^{R}$ would also be smaller. Note that correlations are
higher for the case $L=2$ than for $L=10$. In other words, when the distance
between the escape and target regions is larger, the particle would generally
take longer time to diffuse in the confining domain before the escape and thus
to decorrelate these random variables. Even though all moments vanish as
$r_{0}\to L$, the correlation coefficient gets a nontrivial limit. Curiously,
both curves exhibit slightly non-monotonous behavior with respect to $r_{0}$.

Figure 4:  The correlation coefficient between the escape time $T_{0}$ and the
boundary local time $\ell_{T_{0}}^{R}$ as a function of $r_{0}/L$, for
diffusion between two concentric spheres of radii $R$ and $L$, with $R=1$,
$D=1$, and two values of $L$ indicated in the legend.

###  III.2 First-crossing times

According to Eq. (51), we have

$\tilde{J}_{R}(\ell,p|\bm{x}_{0})=g_{0}(r_{0})\biggl{(}\delta(\ell)-\mu_{0}^{(p)}e^{-\mu_{0}^{(p)}\ell}\biggr{)}$
(82)

and thus

$\tilde{U}(\ell,p|\bm{x}_{0})=g_{0}(r_{0})\,e^{-\mu_{0}^{(p)}\ell}.$ (83)

Figure 5:  (a) Probability flux density $J_{R}(\ell,t|\bm{x}_{0})$ for
diffusion between two concentric spheres of radii $R$ and $L$, with $R=1$,
$L=2$, $r_{0}=1.5$, and $D=1$. Surface shows the regular part of this density,
which was obtained by the numerical Laplace transform inversion of Eq. (82)
with respect to $p$ by Talbot algorithm. Black solid line presents the
prefactor $J_{\infty}^{R}(t|\bm{x}_{0})$ in front of the singular term
$\delta(\ell)$. Note that this curve should be located at $\ell=0$, which is
not visible on the logarithmic scale; it was thus artificially put at
$\ell=10^{-2}$ for illustration purposes. (b) Probability density
$U(\ell,t|\bm{x}_{0})$ of the first-crossing time $\mathcal{T}_{\ell}$ for the
same setting. This density is obtained by the numerical Laplace transform
inversion of Eq. (83) with respect to $p$ by Talbot algorithm. According to
Eq. (55), $J_{R}(\ell,t|\bm{x}_{0})$ is the derivative of
$U(\ell,t|\bm{x}_{0})$ with respect to $\ell$.

Figure 5(a) illustrates the behavior of $J_{R}(\ell,t|\bm{x}_{0})$. As
previously, we plot the regular and singular parts by a surface and a black
curve. In sharp contrast to $J_{D}(\ell,t|\bm{x}_{0})$ shown on Fig. 2, the
regular part of $J_{R}(\ell,t|\bm{x}_{0})$ is negative, as discussed in Sec.
II.5. Moreover, the minimum of the regular part is not shifted with respect to
the maximum of the singular term $J_{\infty}^{R}(t|\bm{x}_{0})$. This is
related to the fact that both extrema are determined by the time needed for
the particle to travel from $\bm{x}_{0}$ to the target region at $R$.
Interestingly, the “amplitude” of the regular part of
$J_{R}(\ell,t|\bm{x}_{0})$ is an order of magnitude higher than that of
$J_{D}(\ell,t|\bm{x}_{0})$. In particular, the black curve showing
$J_{\infty}^{R}(t|\bm{x}_{0})$ looks almost flat at this scale. This
observation does not contradict Eq. (55) that ensures the positivity of the
integral of $J_{R}(\ell^{\prime},t|\bm{x}_{0})$ over $\ell^{\prime}$ from $0$
to any $\ell$. This is confirmed by Fig. 5(b) showing this integral. For each
value of the threshold $\ell$, this figure gives the probability density
$U(\ell,t|\bm{x}_{0})$ of the first-crossing time $\mathcal{T}_{\ell}$. We
recall that this density is not normalized to $1$ due to escape events:

$\displaystyle\int\limits_{0}^{\infty}dt\,U(\ell,t|\bm{x}_{0})$
$\displaystyle=\tilde{U}(\ell,0|\bm{x}_{0})$ (84)
$\displaystyle=\frac{R(L-r_{0})}{r_{0}(L-R)}e^{-\ell/(1/R+1/(L-R))}\leq 1.$

The probability density $U(\ell,t|\bm{x}_{0})$ as a function of both $\ell$
and $t$ exhibits a single “boss”; in fact, it is unlikely to cross a given
threshold $\ell$ at too short or too long times; in turn, the decrease at
large $\ell$ is ensured by the normalization relation (84). Even though the
shape of $U(\ell,t|\bm{x}_{0})$ looks similar to that of
$J_{D}(\ell,t|\bm{x}_{0})$ shown on Fig. 2, their probabilistic
interpretations are different.

Yet another difference between $J_{D}(\ell,t|\bm{x}_{0})$ and
$J_{R}(\ell,t|\bm{x}_{0})$ is that the former strongly depends on the location
of the escape region (compare two panels of Fig. 2), whereas the latter
exhibits only a weak dependence on $L$ whenever $L$ is large enough. For this
reason, we do not present the graphs of $J_{R}(\ell,t|\bm{x}_{0})$ and
$U(\ell,t|\bm{x}_{0})$ for $L=10$, because they are almost indistinguishable
from that shown on Fig. 5 for $L=2$. This can be seen from Eq. (83): when
$\sqrt{p/D}(L-r_{0})\gg 1$, one has

$\tilde{U}(\ell,p|\bm{x}_{0})\approx\frac{R}{r_{0}}e^{-(r_{0}-R+\ell)\sqrt{p/D}-\ell/R}\,,$
(85)

which is independent of $L$. Its inverse Laplace transform yields

$U(\ell,t|\bm{x}_{0})\approx\frac{Re^{-\ell/R}}{r_{0}}\,\frac{(r_{0}-R+\ell)e^{-(r_{0}-R+\ell)^{2}/(4Dt)}}{\sqrt{4\pi
Dt^{3}}}\,.$ (86)

The right-hand side is actually the exact form of the probability density
$U(\ell,t|\bm{x}_{0})$ for a spherical target in the three-dimensional space
(i.e., in the limit $L\to\infty$), as reported in Grebenkov20b . As the
approximation (85) fails in the limit $p\to 0$ (in which case
$\sqrt{p/D}(L-r_{0})\gg 1$ cannot hold), the approximation (86) fails at long
times. This is expected because if the particle has enough time to explore the
domain $\Omega$, it will unavoidably hit the escape region.

Figure 6 allows one to compare the probability density $U(\ell,t|\bm{x}_{0})$
for a particular value $\ell=1$ with the results of Monte Carlo simulations.
One observes an excellent agreement between the theory and simulations. This
figure also highlights the difficulty in getting the values of
$U(\ell,t|\bm{x}_{0})$ at short and long times by Monte Carlo simulations. We
also plot the approximation (86) in the limit $L\to\infty$. As said earlier,
this approximation is very accurate at short times but fails at long times. As
$L$ increases, the validity range of the approximation progressively extends
to longer and longer times.

Figure 6:  Probability density $U(\ell,t|\bm{x}_{0})$ of the first-crossing
time $\mathcal{T}_{\ell}$ of a threshold $\ell$ by the boundary local time
$\ell_{t}^{R}$ for diffusion between two concentric spheres of radii $R$ and
$L$, with $R=1$, $L=2$, $r_{0}=1.5$, $D=1$, and $\ell=1$. Symbols present a
renormalized histogram from Monte Carlo simulations with $M=10^{4}$,
$\delta=10^{-4}$ and $a=5\sigma$, solid line shows the numerical Laplace
transform inversion of Eq. (83) by Talbot algorithm, while dashed line
indicates the approximation (86) corresponding to the limit $L\to\infty$.

##  IV Conclusion

In this paper, we revised the encounter-based approach to diffusion-mediated
surface phenomena and generalized this formalism by allowing a generic
partition of the boundary $\partial\Omega$ into three parts: a target region
$\partial\Omega_{R}$, a reflecting region $\partial\Omega_{N}$ and an escape
region $\partial\Omega_{D}$. While the original formulation in Grebenkov20
dealt with the whole boundary as the target
($\partial\Omega=\partial\Omega_{R}$), our generalization brings a greater
flexibility to modeling of diffusion-controlled reactions and covers a broad
scope of related first-passage problems. From the mathematical point of view,
this extension essentially consists in adding for most involved quantities
Neumann and Dirichlet boundary conditions on $\partial\Omega_{N}$ and
$\partial\Omega_{D}$, respectively. Despite this apparent simplicity, the
inclusion of a “killing mechanism” for the diffusing particle (an escape
through $\partial\Omega_{D}$) has required some conceptual refinements and
appropriate modifications in the formalism. In particular, we focused in the
paper on the probability flux density $j(\bm{x},\ell,t|\bm{x}_{0})$, which was
mostly ignored in former works.

On one hand, the restriction of $j(\bm{x},\ell,t|\bm{x}_{0})$ to the escape
region $\partial\Omega_{D}$ determined, for the first time, the joint
probability density of the escape time $T_{0}$, the particle position
$\bm{X}_{T_{0}}$, and the acquired boundary local time $\ell_{T_{0}}^{R}$ on
the target region $\partial\Omega_{R}$. In particular, the marginal density
$J_{D}(\ell,t|\bm{x}_{0})$, in which the position $\bm{X}_{T_{0}}$ was
averaged out, allows one to characterize not only the statistics of the number
of encounters between the target region and the diffusing particle before its
escape, but also correlations with the escape time. In particular, we obtained
three equivalent representations (44, 48, 49) of the joint moments of
$\ell_{T_{0}}^{R}$ and $T_{0}$ and revealed a probabilistic interpretation of
multiple derivatives of the Green’s function
$\tilde{G}_{q}(\bm{x},p|\bm{x}_{0})$ with respect to $p$ and $q$. Note that
earlier works dealt exclusively with the escape time itself.

On the other hand, the restriction of $j(\bm{x},\ell,t|\bm{x}_{0})$ to the
target region $\partial\Omega_{R}$ required a subtle probabilistic
interpretation in terms of the first-crossing time $\mathcal{T}_{\ell}$ of a
given threshold $\ell$ by the boundary local time $\ell_{t}^{R}$. In
particular, the integral of $j(\bm{x},\ell^{\prime},t|\bm{x}_{0})$ over
$\bm{x}\in\partial\Omega_{R}$ and over $\ell^{\prime}$ from $0$ to $\ell$
determined the probability density $U(\ell,t|\bm{x}_{0})$ of
$\mathcal{T}_{\ell}$. The latter played a key role in the generalization of
conventional surface reactions, described by the Robin boundary condition, to
more general mechanisms. While this generalization was already presented in
Grebenkov20 , its extension in the presence of escape events required some
refinements. For instance, the former derivation relied on the nondecreasing
character of the boundary local time $\ell_{t}$ and its immediate consequence
that $\mathbb{P}\\{\mathcal{T}_{\ell}>t\\}=\mathbb{P}\\{\ell_{t}<\ell\\}$. In
the presence of the escape region, this relation has to be replaced by Eq.
(54) accounting for the value $\mathcal{T}_{\ell}=\infty$ in the case if the
particle has escaped before crossing the threshold $\ell$. More generally, the
lack of a proper normalization of probability densities had to be carefully
managed.

For illustrative purposes, we considered diffusion between concentric spheres,
the inner sphere of radius $R$ being the target region and the outer sphere of
radius $L$ being the escape region. This basic example allowed us to compute
all the discussed quantities in a simple compact form in the Laplace domain
and then to evaluate the inverse Laplace transform numerically. We presented
both $J_{D}(\ell,t|\bm{x}_{0})$ and $J_{R}(\ell,t|\bm{x}_{0})$, as well as
some derived quantities such as $\rho_{D}(\ell|\bm{x}_{0})$ and
$U(\ell,t|\bm{x}_{0})$. We discussed two cases: $L=2$ (when the target and
escape regions are relatively close to each other), and $L=10$ (when they are
well separated). Expectedly, the joint probability density
$J_{D}(\ell,t|\bm{x}_{0})$ strongly depends on the location of the escape
region, whereas $J_{R}(\ell,t|\bm{x}_{0})$ showed only a weak dependence on
$L$. In particular, the probability density $U(\ell,t|\bm{x}_{0})$ approaches
the well-known explicit form (86) for a spherical target in the three-
dimensional space.

While we considered ordinary diffusion inside the confining domain, the
presented formalism allows one to easily incorporate a first-order reaction
kinetics in the bulk. This kinetics may account for a spontaneous
disappearance of the diffusing particle (or its “activity”) with a constant
rate $\gamma$ due to radioactive disintegration, photobleaching, relaxation of
its excited state, disassembly, failure, or biological death. Whatever the
actual killing mechanism, such a “mortal” particle can be treated as having a
random lifetime $\delta$ that obeys the exponential law
$\mathbb{P}\\{\delta>t\\}=e^{-\gamma t}$ Yuste13 ; Meerson15 ; Grebenkov17f ;
Meerson19 . In the conventional approach, the killing mechanism in the bulk is
included by adding the term $-\gamma G_{q}(\bm{x},t|\bm{x}_{0})$ to the right-
hand side of the diffusion equation (1). Its Laplace transform with respect to
$t$ yields the same modified Helmholtz equation (19), in which $p$ is replaced
by $p^{\prime}=p+\gamma$. As a consequence, most results that we obtained in
the Laplace domain, remain valid up to this trivial change. Moreover, many
Laplace-transformed quantities, evaluated at $p=\gamma$, admit useful
probabilistic interpretations in terms of the stopping condition at the
“death” time $\delta$. For instance,

$\gamma\,\tilde{\rho}(\ell,\gamma|\bm{x}_{0})=\int\limits_{0}^{\infty}dt\,\gamma
e^{-\gamma t}\,\rho(\ell,t|\bm{x}_{0})$ (87)

is the probability density of the boundary local time $\ell_{\delta}^{R}$ at
the random time $\delta$ of the particle “death” in the bulk, with $\gamma
e^{-\gamma t}$ being the probability density of $\delta$. In other words, many
Laplace-transformed quantities that we derived in this paper, have their own
interest, even without evaluating their inverse Laplace transforms. The
combined effect of the first-order kinetics in the bulk, surface reactions on
the target region, and escape events can be further explored.

The present work can be extended in several directions. On the mathematical
side, the spectral properties of the Dirichlet-to-Neumann operator need
further attention. In addition to a rigorous demonstration of the announced
basic properties, the asymptotic behavior of the eigenvalues $\mu_{k}^{(p)}$
has to be uncovered, in particular, in the limit when either
$\partial\Omega_{R}$ or $\partial\Omega_{D}$ (or both) is small. We expect
that some matched asymptotic tools Bressloff22a ; Bressloff22b  can be adapted
to investigate this problem. On the application side, one can investigate in
more detail the effect of escape events onto various surface reaction
mechanisms introduced in Grebenkov20 . The proposed extension can be further
analyzed in the presence of multiple independently diffusing particles
Grebenkov22b , eventual resetting mechanisms Evans20 ; Bressloff22c ;
Benkhadaj22  and drifts Grebenkov22a . Future applications of the extended
encounter-based approach can bring more realistic features to a theoretical
description of biology-inspired transport processes, notably in living cells.

###### Acknowledgements.

The author thanks the Alexander von Humboldt Foundation for support through
the Bessel Prize award.

## References

  * (1) B. Alberts, A. Johnson, J. Lewis, D. Morgan, M. Raff, K. Roberts, and P. Walter, Molecular Biology of the Cell (Garland Science, New York, NY, 2014).
  * (2) S. A. Rice, Diffusion-limited reactions (Elsevier, Amsterdam, 1985).
  * (3) D. A. Lauffenburger and J. Linderman, Receptors: Models for Binding, Trafficking, and Signaling (Oxford University Press, Oxford, 1993).
  * (4) Z. Schuss, Brownian Dynamics at Boundaries and Interfaces in Physics, Chemistry and Biology (Springer, New York, 2013).
  * (5) K. Lindenberg, R. Metzler, and G. Oshanin Chemical Kinetics: Beyond the Textbook (World Scientific, New Jersey, 2019).
  * (6) P. C. Bressloff and J. M. Newby, Stochastic models of intracellular transport, Rev. Mod. Phys. 85, 135-196 (2013).
  * (7) S. Redner, A Guide to First Passage Processes (Cambridge, Cambridge University press, 2001).
  * (8) R. Metzler, G. Oshanin, and S. Redner (Eds), First-Passage Phenomena and Their Applications (Singapore, World Scientific, 2014).
  * (9) S. Condamin, O. Bénichou, V. Tejedor, R. Voituriez, and J. Klafter, First-passage time in complex scale-invariant media, Nature 450, 77 (2007).
  * (10) O. Bénichou, C. Chevalier, J. Klafter, B. Meyer, and R. Voituriez, Geometry-controlled kinetics, Nat. Chem. 2, 472-477 (2010).
  * (11) D. Holcman and Z. Schuss, Control of flux by narrow passages and hidden targets in cellular biology, Phys. Progr. Rep. 76, 074601 (2013).
  * (12) O. Bénichou and R. Voituriez, From first-passage times of random walks in confinement to geometry-controlled kinetics, Phys. Rep. 539, 225-284 (2014).
  * (13) D. S. Grebenkov, Universal formula for the mean first passage time in planar domains, Phys. Rev. Lett. 117, 260201 (2016).
  * (14) T. Guérin, N. Levernier, O. Bénichou, and R. Voituriez, Mean first-passage times of non-Markovian random walkers in confinement, Nature 534, 356-359 (2016).
  * (15) Y. Lanoiselée, N. Moutal, and D. S. Grebenkov, Diffusion-limited reactions in dynamic heterogeneous media, Nature Commun. 9, 4398 (2018).
  * (16) C. W. Gardiner, Handbook of stochastic methods for physics, chemistry and the natural sciences (Springer: Berlin, 1985).
  * (17) N. G. Van Kampen, Stochastic Processes in Physics and Chemistry (Elsevier, Amsterdam, 1992).
  * (18) F. C. Collins and G. E. Kimball, Diffusion-controlled reaction rates, J. Colloid Sci. 4, 425 (1949).
  * (19) H. C. Berg and E. M. Purcell, Physics of chemoreception, Biophys. J. 20, 193 (1977).
  * (20) H. Sano and M. Tachiya, Partially diffusion-controlled recombination, J. Chem. Phys. 71, 1276 (1979).
  * (21) K. R. Brownstein and C. E. Tarr, Importance of Classical Diffusion in NMR Studies of Water in Biological Cells, Phys. Rev. A 19, 2446-2453 (1979).
  * (22) G. H. Weiss, Overview of theoretical models for reaction rates, J. Stat. Phys. 42, 3 (1986).
  * (23) J. G. Powles, M. J. D. Mallett, G. Rickayzen, and W. A. B. Evans, Exact analytic solutions for diffusion impeded by an infinite array of partially permeable barriers, Proc. R. Soc. London A 436, 391-403 (1992).
  * (24) B. Sapoval, General Formulation of Laplacian Transfer Across Irregular Surfaces, Phys. Rev. Lett. 73, 3314 (1994).
  * (25) B. Sapoval, M. Filoche, and E. Weibel, Smaller is better – but not too small: A physical scale for the design of the mammalian pulmonary acinus, Proc. Nat. Ac. Sci. USA 99, 10411-10416 (2002).
  * (26) D. S. Grebenkov, M. Filoche, B. Sapoval, and M. Felici, Diffusion-reaction in Branched Structures: Theory and Application to the Lung Acinus, Phys. Rev. Lett. 94, 050602 (2005).
  * (27) S. D. Traytak and W. Price, Exact solution for anisotropic diffusion-controlled reactions with partially reflecting conditions, J. Chem. Phys. 127, 184508 (2007).
  * (28) P. C. Bressloff, B. A. Earnshaw, and M. J. Ward, Diffusion of protein receptors on a cylindrical dendritic membrane with partially absorbing traps, SIAM J. Appl. Math. 68, 1223-1246 (2008).
  * (29) A. Singer, Z. Schuss, A. Osipov, and D. Holcman, Partially reflected diffusion, SIAM J. Appl. Math. 68, 844 (2008).
  * (30) D. S. Grebenkov, Searching for partially reactive sites: Analytical results for spherical targets, J. Chem. Phys. 132, 034104 (2010).
  * (31) S. D. Lawley and J. P. Keener, A new derivation of Robin boundary conditions through homogenization of a stochastically switching boundary, SIAM J. Appl. Dyn. Syst. 14, 1845-1867 (2015).
  * (32) D. S. Grebenkov, Analytical representations of the spread harmonic measure density, Phys. Rev. E 91, 052108 (2015).
  * (33) A. S. Serov, C. Salafia, D. S. Grebenkov, and M. Filoche, The Role of Morphology in Mathematical Models of Placental Gas Exchange, J. Appl. Physiol. 120, 17-28 (2016).
  * (34) P. C. Bressloff, Stochastic switching in biology: from genotype to phenotype, J. Phys. A. 50, 133001 (2017).
  * (35) D. S. Grebenkov and G. Oshanin, Diffusive escape through a narrow opening: new insights into a classic problem, Phys. Chem. Chem. Phys. 19, 2723-2739 (2017).
  * (36) F. Piazza and D. S. Grebenkov, Diffusion-controlled reaction rate on non-spherical partially absorbing axisymmetric surfaces, Phys. Chem. Chem. Phys. 21, 25896 (2019).
  * (37) D. S. Grebenkov, Spectral theory of imperfect diffusion-controlled reactions on heterogeneous catalytic surfaces, J. Chem. Phys. 151, 104108 (2019).
  * (38) H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids, 2nd Ed. (Oxford University Press, 1959).
  * (39) J. Crank, The Mathematics of Diffusion (Oxford University Press, 1956).
  * (40) R. K. M. Thambynayagam, The Diffusion Handbook: Applied Solutions for Engineers (McGraw Hill, 2011).
  * (41) D. S. Grebenkov and B.-T. Nguyen, Geometrical structure of Laplacian eigenfunctions, SIAM Rev. 55, 601-667 (2013).
  * (42) V. G. Maz’ya, S. A. Nazarov, and B. A. Plamenevskii, Asymptotic expansions of the eigenvalues of boundary value problems for the Laplace operator in domains with small holes, Math. USSR. Izv 24, 321-345 (1985).
  * (43) M. J. Ward and J. B. Keller, Strong localized perturbations of eigenvalue problems, SIAM J. Appl. Math. 53, 770-798 (1993).
  * (44) T. Kolokolnikov, M. S. Titcombe, and M. J. Ward, Optimizing the fundamental Neumann eigenvalue for the Laplacian in a domain with small traps, Eur. J. Appl. Math 16, 161 (2005).
  * (45) A. Singer, Z. Schuss, D. Holcman, and R. S. Eisenberg, Narrow escape, part I, J. Stat. Phys. 122, 437-463 (2006).
  * (46) A. Singer, Z. Schuss, and D. Holcman, Narrow escape, part II: the circular disk, J. Stat. Phys. 122, 465 (2006).
  * (47) A. Singer, Z. Schuss, and D. Holcman, Narrow escape, part III: non-smooth domains and Riemann surfaces, J. Stat. Phys. 122, 491 (2006).
  * (48) Z. Schuss, A. Singer, and D. Holcman, The narrow escape problem for diffusion in cellular microdomains, Proc. Nat. Acad. Sci. USA 104, 16098-16103 (2007).
  * (49) O. Bénichou and R. Voituriez, Narrow-escape time problem: time needed for a particle to exit a confining domain through a small window, Phys. Rev. Lett. 100, 168105 (2008).
  * (50) S. Pillay, M. J. Ward, A. Peirce, and T. Kolokolnikov, An asymptotic analysis of the mean first passage time for narrow escape problems: part I: two-dimensional domains, Multiscale Model. Simul. 8, 803-835 (2010).
  * (51) A. F. Cheviakov, M. J. Ward, and R. Straube, An asymptotic analysis of the mean first passage time for narrow escape problems: part II: the sphere, Multiscale Model. Simul. 8, 836-870 (2010).
  * (52) A. F. Cheviakov and M. J. Ward, Optimizing the principal eigenvalue of the Laplacian in a sphere with interior traps, Math. Comput. Modelling 53, 1394-1409 (2011).
  * (53) A. F. Cheviakov, A. S. Reimer, and M. J. Ward, Mathematical modeling and numerical computation of narrow escape problems, Phys. Rev. E 85, 021131 (2012).
  * (54) D. Holcman and Z. Schuss, The narrow escape problem, SIAM Rev. 56, 213-257 (2014).
  * (55) T. Agranov and B. Meerson, Narrow escape of interacting diffusing particles, Phys. Rev. Lett. 120, 120601 (2018).
  * (56) D. S. Grebenkov, R. Metzler, and G. Oshanin, Full distribution of first exit times in the narrow escape problem, New J. Phys. 21, 122001 (2019).
  * (57) S. D. Traytak, The diffusive interaction in diffusion-limited reactions: the steady-state case, Chem. Phys. Lett. 197, 247-254 (1992).
  * (58) S. Condamin and O. Bénichou, Exact expressions of mean first-passage times and splitting probabilities for random walks in bounded rectangular domains, J. Chem. Phys. 124, 206103 (2006).
  * (59) C. Chevalier, O. Bénichou, B. Meyer and R. Voituriez, First-passage quantities of Brownian motion in a bounded domain with multiple targets: a unified approach, J. Phys. A: Math. Theor. 44, 025002 (2011).
  * (60) M. Galanti, D. Fanelli, S. D. Traytak, and F. Piazza, Theory of diffusion-influenced reactions in complex geometries Phys. Chem. Chem. Phys. 18, 15950-15954 (2016).
  * (61) D. S. Grebenkov and S. Traytak, Semi-analytical computation of Laplacian Green functions in three-dimensional domains with disconnected spherical boundaries, J. Comput. Phys. 379, 91-117 (2019).
  * (62) D. S. Grebenkov, Diffusion toward non-overlapping partially reactive spherical traps: fresh insights onto classic problems, J. Chem. Phys. 152, 244108 (2020).
  * (63) J. Klinger, R. Voituriez, and O. Bénichou, Splitting Probabilities of Symmetric Jump Processes, Phys. Rev. Lett. 129, 140603 (2022).
  * (64) D. S. Grebenkov, Paradigm Shift in Diffusion-Mediated Surface Phenomena, Phys. Rev. Lett. 125, 078102 (2020).
  * (65) P. Lévy, Processus Stochastiques et Mouvement Brownien (Paris: Gauthier-Villard, 1965).
  * (66) K. Itô and H. P. McKean, Diffusion Processes and Their Sample Paths (Berlin: Springer, 1965).
  * (67) M. Freidlin, Functional Integration and Partial Differential Equations (Annals of Mathematics Studies, Princeton, NJ: Princeton University Press, 1985).
  * (68) A. N. Borodin and P. Salminen, Handbook of Brownian Motion: Facts and Formulae (Birkhauser Verlag, Basel-Boston-Berlin, 1996).
  * (69) S. N. Majumdar, Brownian functionals in physics and computer science, Curr. Sci. 88, 2076-2092 (2005).
  * (70) D. S. Grebenkov, Statistics of boundary encounters by a particle diffusing outside a compact planar domain, J. Phys. A.: Math. Theor. 54, 015003 (2021).
  * (71) D. S. Grebenkov, Joint distribution of multiple boundary local times and related first-passage time problems with multiple targets, J. Stat. Mech. 103205 (2020).
  * (72) W. Arendt, A. F. M. ter Elst, J. B. Kennedy, and M. Sauter, The Dirichlet-to-Neumann operator via hidden compactness, J. Funct. Anal. 266, 1757-1786 (2014).
  * (73) D. Daners, Non-positivity of the semigroup generated by the Dirichlet-to-Neumann operator, Positivity 18, 235-256 (2014).
  * (74) A. F. M. ter Elst and E. M. Ouhabaz, Analysis of the heat kernel of the Dirichlet-to-Neumann operator, J. Funct. Anal. 267, 4066-4109 (2014).
  * (75) J. Behrndt and A. F. M. ter Elst, Dirichlet-to-Neumann maps on bounded Lipschitz domains, J. Differ. Equ. 259, 5903-5926 (2015).
  * (76) W. Arendt and A. F. M. ter Elst, The Dirichlet-to-Neumann operator on exterior domains, Potential Anal. 43, 313-340 (2015).
  * (77) A. Hassell and V. Ivrii, Spectral asymptotics for the semiclassical Dirichlet to Neumann operator J. Spectr. Theory 7, 881-905 (2017).
  * (78) A. Girouard and I. Polterovich, Spectral geometry of the Steklov problem, J. Spectr. Theory 7, 321-359 (2017).
  * (79) D. S. Grebenkov, Probability distribution of the boundary local time of reflected Brownian motion in Euclidean domains, Phys. Rev. E 100, 062110 (2019).
  * (80) P. Mörters and Y. Peres, Brownian Motion (Cambridge Series in Statistical and Probabilistic Mathematics, Cambridge University Press, New York, 2010).
  * (81) D. S. Grebenkov, Surface Hopping Propagator: An Alternative Approach to Diffusion-Influenced Reactions, Phys. Rev. E 102, 032125 (2020).
  * (82) D. S. Grebenkov, Statistics of diffusive encounters with a small target: Three complementary approaches, J. Stat. Mech. 083205 (2022).
  * (83) A. Talbot, The accurate numerical inversion of Laplace transforms, J. Inst. Math. Appl. 23, 97-120 (1979).
  * (84) A. Godec and R. Metzler, Universal Proximity Effect in Target Search Kinetics in the Few-Encounter Limit, Phys. Rev. X 6, 041037 (2016).
  * (85) D. S. Grebenkov, R. Metzler, and G. Oshanin, Strong defocusing of molecular reaction times results from an interplay of geometry and reaction control, Commun. Chem. 1, 96 (2018).
  * (86) S. B. Yuste, E. Abad, and K. Lindenberg, Exploration and trapping of mortal random walkers, Phys. Rev. Lett. 110, 220603 (2013).
  * (87) B. Meerson and S. Redner, Mortality, Redundancy, and Diversity in Stochastic Search, Phys. Rev. Lett. 114, 198101 (2015).
  * (88) D. S. Grebenkov and J.-F. Rupprecht, The escape problem for mortal walkers, J. Chem. Phys. 146, 084106 (2017).
  * (89) B. Meerson, Mortal Brownian motion: Three short stories, Int. J. Mod. Phys. B 33, 1950172 (2019).
  * (90) P. C. Bressloff, Diffusion-mediated absorption by partially-reactive targets: Brownian functionals and generalized propagators, J. Phys. A: Math. Theor. 55, 205001 (2022).
  * (91) P. C. Bressloff, Narrow capture problem: an encounter-based approach to partially reactive targets, Phys. Rev. E 105, 034141 (2022).
  * (92) D. S. Grebenkov, Depletion of Resources by a Population of Diffusing Species, Phys. Rev. E 105, 054402 (2022).
  * (93) M. R. Evans, S. N. Majumdar, and G. Schehr, Stochastic resetting and applications J. Phys. A: Math. Theor. 53, 193001 (2020).
  * (94) P. C. Bressloff, Diffusion-mediated surface reactions and stochastic resetting, J. Phys. A: Math. Theor. 55, 275002 (2022).
  * (95) Z. Benkhadaj and D. S. Grebenkov, Encounter-based approach to diffusion with resetting, Phys. Rev. E 106, 044121 (2022).
  * (96) D. S. Grebenkov, An encounter-based approach for restricted diffusion with a gradient drift, J. Phys. A: Math. Theor. 55, 045203 (2022).