# Escape-from-a-layer approach for simulating the boundary local time in Euclidean domains

Yilin Ye, Adrien Chaigneau, Denis S Grebenkov

► To cite this version:

Yilin Ye, Adrien Chaigneau, Denis S Grebenkov. Escape-from-a-layer approach for simulating the boundary local time in Euclidean domains. Journal of Computational Physics, 2025, 537, pp.114099. 10.1016/j.jcp.2025.114099. hal-05344181

HAL Id: hal-05344181

https://hal.science/hal-05344181v1

Submitted on 3 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Escape-from-a-layer approach for simulating the boundary local time in Euclidean domains

Yilin Ye<sup>a</sup>, Adrien Chaigneau<sup>a</sup>, Denis S. Grebenkov<sup>a</sup>

<sup>a</sup>Laboratoire de Physique de la Matière Condensée,
CNRS – Ecole Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

**Abstract**

We propose an efficient numerical approach to simulate the boundary local time of reflected Brownian motion, as well as the time and position of the associated reaction event on a smooth boundary of a Euclidean domain. This approach combines the standard walk-on-spheres algorithm in the bulk with the approximate solution of the escape problem in a boundary layer. In this way, the most time-consuming simulation of multiple reflections on the boundary is replaced by an equivalent escape event. We validate the proposed escape-from-a-layer approach by comparing simulated statistics of the boundary local time with exact results known for simple domains (a disk, a circular annulus, a sphere, a spherical shell) and with the numerical results obtained by a finite-element method in more sophisticated domains. This approach offers a powerful tool for simulating reflected Brownian motion in multi-scale confinements such as porous media or biological environments, and for solving the related partial differential equations. Its applications in the context of diffusion-controlled reactions in chemical physics are discussed.

**Keywords:** boundary local time, reflected Brownian motion, escape problem, walk-on-spheres, Monte Carlo simulation, diffusion-controlled reaction, Robin and Neumann boundary conditions

**1. Introduction**

The boundary local time $\ell_t$ plays the central role in the theory of stochastic processes [1–3]. For instance, reflected Brownian motion $X_t$ inside a given Euclidean domain $\Omega \subset \mathbb{R}^d$ with a smooth boundary $\partial\Omega$ can be constructed as the solution of the Skorokhod stochastic equation [2, 3]

$$dX_t = \sqrt{2D} dW_t + n(X_t) d\ell_t, \quad X_0 = x_0, \quad \ell_0 = 0, \quad (1)$$

where $W_t$ is the standard Wiener process in $\mathbb{R}^d$, $D$ is the diffusion coefficient, $x_0 \in \Omega$ is the starting point, $n(x)$ is the unit normal vector at a boundary point $x \in \partial\Omega$ oriented inward the domain $\Omega$, and $\ell_t$ is a nondecreasing stochastic process that increments at each encounter of $X_t$ with the boundary. Qualitatively, the first term in Eq. (1) describes ordinary Brownian motion inside $\Omega$, whereas the second term ensures that $X_t$ is reflected back normally into $\Omega$ at each encounter with the boundary. From the physical point of view, one can think of the second term as an infinitely local force field that pushes the diffusing particle back to the confining domain [4]. Curiously, the single stochastic equation (1) determines simultaneously two tightly related stochastic processes: the position $X_t$ and the boundary local time $\ell_t$ (which has, despite its name, units of length). The latter should not be confused with a point local time, which represents the residence time of Brownian motion in a vicinity of a bulk point (see [5–7]). Note that the boundary local time can be expressed as

$$\ell_t = \lim_{\epsilon \to 0} \frac{D}{\epsilon} \int_0^t dt' \Theta(\epsilon - |X_{t'} - \partial\Omega|), \quad (2)$$

where $|x - \partial\Omega|$ is the Euclidean distance between a point $x$ and the boundary $\partial\Omega$, and $\Theta(z)$ is the Heaviside step function: $\Theta(z) = 1$ for $z > 0$ and 0 otherwise. The integral represents the residence time of the process $X_t$ in a thin
layer $\partial\Omega_{\epsilon} = \{x \in \Omega : |x - \partial\Omega| < \epsilon\}$ of width $\epsilon$ near the boundary $\partial\Omega$, which is rescaled by $\epsilon$ to get a nontrivial limit, as $\epsilon \to 0$. Alternatively, one has

$$ \ell_t = \lim_{\epsilon \to 0} \epsilon N_t^{(\epsilon)}, \qquad (3) $$

where $N_t^{(\epsilon)}$ is the number of crossings of the layer $\partial\Omega_{\epsilon}$ up to time $t$. In other words, the boundary local time $\ell_t$ can be understood either as the rescaled residence time in a thin boundary layer, or as the rescaled number of encounters with that boundary.

Apart from its major role in the theory of stochastic processes, the boundary local time was recently employed as the conceptual pillar to build the encounter-based approach to diffusion-controlled reactions [4]. This approach laid a theoretical ground for assessing the statistics of encounters between the diffusing particle and the boundary, provided an intuitively clear probabilistic interpretation of partial reactivity and allowed one to consider much more general surface reaction mechanisms [8–15]. In fact, the reaction event of a diffusing particle on the boundary can be initiated when the boundary local time $\ell_t$ exceeds a chosen random threshold $\hat{\ell}$. In this way, the probability distribution of the threshold determines the surface reaction mechanism *independently* of the diffusive dynamics. So, the exponential distribution of $\hat{\ell}$ corresponds to the standard constant reactivity implemented via Robin boundary condition, whereas other distributions can describe, for instance, progressive activation or passivation of the boundary by the diffusing particle [15, 16]. While the former setting can be simulated by conventional methods, implementation of more sophisticated surface reaction mechanisms requires the knowledge of the boundary local time. Moreover, these concepts were further extended to describe permeation across membranes [17–19], diffusive exchange between adjacent compartments [20], and the related snapping out Brownian motion [21, 22].

Despite these advances, the fundamental relation between the statistics of the boundary local time $\ell_t$ and the geometrical structure of the confining domain $\Omega$ remains poorly understood. For instance, the distribution of $\ell_t$ is known explicitly only for two simple domains: a half-line and the exterior of a sphere [8, 23, 24]. Indeed, the symmetries of these domains allow for the separation of variables and thus reduction to a one-dimensional problem that can be solved exactly. In contrast, most diffusive processes in nature and industrial applications occur in multi-scale media with irregular, sophisticated boundaries, for which numerical tools are needed to access the statistics of encounters. Even though the spectral expansions developed in [4] are formally still applicable in such domains, their numerical implementation (e.g., by a finite-element method) can be very costly and time-consuming, especially in three dimensions. In this light, one may prefer Monte Carlo techniques that offer great flexibility and moderate computational costs [25–27]. While these techniques are broadly used for simulating stochastic processes and diffusion-controlled reactions in complex domains [28–50], their adaptation for studying the boundary local time is not straightforward.

The most basic way to simulate the boundary local time consists in approximating reflected Brownian motion by a random walk on a lattice with spacing $\epsilon$ and counting the number $N_t^{(\epsilon)}$ of its encounters with the discretized boundary $\partial\Omega$. According to Eq. (3), $\epsilon N_t^{(\epsilon)}$ is an approximation of $\ell_t$, if $\epsilon$ is small enough. The need for boundary discretization and excessively long trajectories (with $2dDt/\epsilon^2$ steps) are major drawbacks of this method. The first drawback can be relaxed by performing off-lattice random walks, i.e., a sequence of centered Gaussian jumps of variance $\sigma^2$ with normal reflections on the boundary. In this case, one can either count the number $N_t^{(\sigma)}$ of reflections and use again Eq. (3) with $\sigma$ instead of $\epsilon$, or calculate the residence time of this random walk inside a boundary layer of width $\epsilon \sim \sigma$ and use Eq. (2) to approximate $\ell_t$ (see also [51] for the discussion of Euler schemes). Since a typical one-jump displacement $\sigma$ should be the smallest length scale of the problem, this method can be efficient for simple domains (see, e.g., [52]) but it becomes too time-consuming in multi-scale media. This issue was partly resolved by Zhou et al. who combined the standard walk-on-spheres (WOS) algorithm by Muller [53] and the fixed-length displacements in a thin layer near boundary [54]. Their hybrid method exploits the well-known advantages of the WOS algorithm for fast simulations of large-scale displacements far from the boundary (see details in Sec. 2.1). As a consequence, the major remaining computational limitation is caused by modeling small displacements inside the boundary layer of small width $\epsilon$. This method was further improved and adapted to simulating snapping out Brownian motion in planar bounded domains by Schumm and Bressloff [55]. In particular, they suggested using walk-on-spheres even inside the boundary layer of width $\epsilon$, with the jump distance being fixed to be $2\epsilon$, in order to speed up simulations within that layer. Despite this progress, an accurate estimation of the boundary local time $\ell_t$ within a boundary layer requires taking $\epsilon$ sufficient small so that modeling of reflected Brownian motion within this layer remains the critical bottleneck of these multi-scale Monte Carlo simulations.

In this paper, we propose a different approach to simulate the stochastic process $(X_t, \ell_t)$ inside a confining domain.
<figure><img src="image_2.png" /><figcaption>Figure 1: (a) A zoom near the boundary $\partial\Omega$: once the WOS arrives inside the boundary layer $\partial\Omega_{\varepsilon}$ of width $\varepsilon$, one aims at replacing a costly detailed simulation of the random trajectory (in green) from the entrance point $\mathbf{x}_k$ to the exit point $\mathbf{x}_{k+1}$ by a single escape event. For this purpose, one needs to generate the escape position $\mathbf{x}_{k+1}$, the escape time $\tau$, and the acquired boundary local time $l'_{\tau}$. (b) Flat layer approximation: starting from a point $(0, y_0)$ (yellow circle), a random trajectory (in blue) experiences numerous reflections on the flat boundary before escaping an infinite stripe of width $\varepsilon$ at a random escape time $\tau$ in a random escape position $(x_{\tau}, \varepsilon)$ (red circle). These reflections are characterized by the acquired boundary local time $l'_{\tau}$. The random trajectory is thus replaced by a simple jump from $(0, y_0)$ to $(x_{\tau}, \varepsilon)$.</figcaption></figure>

We still employ the WOS algorithm for fast simulations of large-scale displacements far from the boundary. In turn, we replace the detailed time-consuming simulation of the stochastic process within a boundary layer by a single escape event from that layer. A similar approach was recently implemented for simulating the escape of a sticky particle [56, 57]. In fact, once the process arrives inside the layer, we introduce the random escape time $\tau$ from that layer, i.e., the first-passage time to the set $\Gamma_{\varepsilon} = \{\mathbf{x} \in \Omega : |\mathbf{x} - \partial\Omega| = \varepsilon\}$ of points in $\Omega$ that are equidistant to the boundary (Fig. 1a).

To realize an escape event, the central quantity of interest is the joint probability density of the escape time $\tau$, the escape position $X_{\tau}$, and the acquired boundary local time $l_{\tau}$. Its exact spectral expansion was recently derived by means of the encounter-based approach [58]. As a smooth boundary is locally flat, the spectral expansion takes a simple form that allows its direct implementation in Monte Carlo simulations. Moreover, when the boundary is composed of flat elements, the simulation of the escape event becomes almost exact. In addition, we propose an improvement to account for the curvature of the boundary in order to enable using a larger width $\varepsilon$ of the boundary layer. In this way, there is no need for simulating multiple reflections on the boundary that was the most time-consuming step and the main source of accumulating errors in former techniques.

The paper is organized as follows. In Sec. 2, we describe the escape-from-a-layer (EFL) approach. We start with a simpler approximation of a flat boundary layer and then discuss its improvement to curved boundary layers. Section 3 summarizes an algorithmic implementation of the EFL approach. Section 4 validates this method for several confining domains by comparing simulated results to either exact values, or numerical values obtained by a finite-element method. We also compare the performance of our method to the state-of-the-art method by Schumm and Bressloff [55]. In Sec. 5, we discuss further improvements, extensions, applications and future perspectives of this method.

## 2. Escape-from-a-layer approach

### 2.1. Walk-on-spheres algorithm

We start by recalling the walk-on-spheres algorithm for simulating Brownian motion in a confining domain $\Omega \subset \mathbb{R}^d$ with an *absorbing* boundary $\partial\Omega$ [53]. From a prescribed starting point $\mathbf{x}_0 \in \Omega$, one determines the Euclidean distance $\rho_0 = |\mathbf{x}_0 - \partial\Omega|$ to the boundary $\partial\Omega$ and draws a ball $B_{\rho_0}(\mathbf{x}_0)$ of radius $\rho_0$ centered at $\mathbf{x}_0$. The continuity of Brownian motion implies that any trajectory must cross the boundary of this ball, $\partial B_{\rho_0}(\mathbf{x}_0)$, before hitting the boundary $\partial\Omega$. One can therefore replace lengthy simulations of the trajectory inside the ball by an escape event, i.e., a jump to a uniformly distributed point $\mathbf{x}_1$ on $\partial B_{\rho_0}(\mathbf{x}_0)$. The probability distribution of the time $\tau_1$ needed to escape from the ball is known explicitly (see Appendix A) that allows one to easily generate the random variable $\tau_1$. From the point
$\mathbf{x}_1$, the previous step is repeated, i.e., one determines the distance $\rho_1 = |\mathbf{x}_1 - \partial\Omega|$ and generates the escape time $\tau_2$ and the escape position $\mathbf{x}_2$ as a uniformly distributed point on the sphere $\partial B_{\rho_1}(\mathbf{x}_1)$, and so on. The algorithm is stopped when the distance to the boundary becomes smaller than a prescribed threshold $\varepsilon$. This procedure samples a sequence of points $X_{t_1}, X_{t_2}, \dots, X_{t_k}$ of a random trajectory $X_t$ at random times $t_1 = \tau_1, t_2 = t_1 + \tau_2, \dots, t_k = t_{k-1} + \tau_k$. As the jump length is progressively adapted according to the proximity to the boundary, large jumps can be performed far from the boundary. This variability ensures remarkably fast convergence of the method [59] and allows for efficient simulations of Brownian motion in multi-scale environments [29–32, 35–38].

The threshold $\varepsilon$ determines the width of a boundary layer, which is considered as absorbing: the process is stopped once it arrives to this layer. The value of $\varepsilon$ can be either assigned from the underlying microscopic model (e.g., the interaction range between the diffusing particle and the boundary), or treated as a numerical parameter of the method. As shown by Binder and Braverman [59], the computational time scales linearly with $\ln(1/\varepsilon)$, i.e., one can choose an extremely small $\varepsilon$ with almost no computational cost. This remarkably fast convergence explains the efficiency and popularity of the WOS algorithm in domains with absorbing boundaries. However, if the boundary is not absorbing, one loses such fast convergence. Indeed, if the process is not stopped upon the first arrival into the boundary layer, one needs to simulate reflections on the boundary that forms the critical bottleneck for an extension of the WOS method to reflecting boundaries. In the next subsections, we propose an efficient way to solve this problem.

## 2.2. The layer escape problem

One can run the WOS algorithm until some time $t_k$ when the process arrives to a point $\mathbf{x}_k$ inside the boundary layer $\Omega_\varepsilon$ of width $\varepsilon$: $\Omega_\varepsilon = \{\mathbf{x} \in \Omega : |\mathbf{x} - \partial\Omega| < \varepsilon\}$. Our goal is to substitute a lengthy simulation of the trajectory $X'_t = X_{t+t_k}$ inside this layer by a single displacement out of $\Omega_\varepsilon$ (as reflected Brownian motion is Markovian, the trajectory $X_t$ with $t < t_k$ does not matter so that we can consider $t_k$ as a new origin of time for the shifted in time process $X'_t$, denoted by prime). For this purpose, we aim at generating three random variables: (i) the escape time $\tau = \inf\{t > 0 : X'_t \notin \Omega_\varepsilon | X'_0 = \mathbf{x}_k\}$ from the layer, i.e., the first instance when the process leaves $\Omega_\varepsilon$ by crossing the surface $\Gamma_\varepsilon = \{\mathbf{x} \in \Omega : |\mathbf{x} - \partial\Omega| = \varepsilon\}$; (ii) the position $X'_\tau$ at the escape time, and (iii) the boundary local time $\ell'_\tau$ acquired during this escape event (Fig. 1a). If we manage to generate these random variables, one can execute a single escape event from the layer, without simulating the trajectory inside it. After the escape, the new position is $X_{t_{k+1}} = \mathbf{x}_{k+1} = X'_\tau$, the time is incremented by $\tau$, $t_{k+1} = t_k + \tau$, and the boundary local time is incremented by $\ell'_\tau$: $\ell_{t_{k+1}} = \ell_{t_k} + \ell'_\tau$. From this position, the WOS algorithm is resumed, until the particle arrives again inside the layer, and so on. The simulation is stopped according to a prescribed stopping condition (see Sec. 2.5).

This layer escape problem was solved in [58]. In fact, for any bounded domain with a smooth boundary, a formal spectral expansion was derived for the probability density flux $j(\mathbf{x}, \ell, t|\mathbf{x}_k)$, which is the joint probability density of the escape position $X'_\tau \in \Gamma_\varepsilon$, the acquired boundary local time $\ell'_\tau$, and the escape time $\tau$. In the Laplace domain, this spectral expansion reads for any $\mathbf{x} \in \Gamma_\varepsilon$:

$$ 
\begin{aligned}
\tilde{j}(\mathbf{x}, \ell, p|\mathbf{x}_k) &= \int_0^\infty dt e^{-pt} j(\mathbf{x}, \ell, t|\mathbf{x}_k) \\
&= \tilde{j}_0(\mathbf{x}, p|\mathbf{x}_k)\delta(\ell) - \sum_{k=0}^\infty [V_k^{(p)}(\mathbf{x}_k)]^* (\partial_n V_k^{(p)}(\mathbf{x}))e^{-\mu_k^{(p)}\ell},
\end{aligned}
\quad (4) $$

where $\partial_n$ is the normal derivative on the boundary $\Gamma_\varepsilon$ oriented outwards the domain $\Omega_\varepsilon$, asterisk denotes complex conjugate, $p \ge 0$ is a scalar, $\tilde{j}_0(\mathbf{x}, p|\mathbf{x}_k)$ is the Laplace transform of the probability density flux onto $\Gamma_\varepsilon$ in the case of the absorbing boundary $\partial\Omega$, while $\mu_k^{(p)}$ and $V_k^{(p)}(\mathbf{x})$ with $k = 0, 1, 2, \dots$ are the eigenvalues and eigenfunctions of the generalized Steklov problem in the boundary layer $\Omega_\varepsilon$:

$$ (p - D\Delta)V_k^{(p)}(\mathbf{x}) = 0 \quad (\mathbf{x} \in \Omega_\varepsilon), \quad (5a) $$

$$ \partial_n V_k^{(p)}(\mathbf{x}) = \mu_k^{(p)} V_k^{(p)}(\mathbf{x}) \quad (\mathbf{x} \in \partial\Omega), \quad (5b) $$

$$ V_k^{(p)}(\mathbf{x}) = 0 \quad (\mathbf{x} \in \Gamma_\varepsilon) \quad (5c) $$

(for mathematical details about the Steklov problem, see [60] and references therein). The first term in Eq. (4) represents the contribution of trajectories that do not hit the boundary $\partial\Omega$ until escaping the layer $\Omega_\varepsilon$. As a consequence,
the acquired boundary local time $\ell'_\tau$ is zero, as represented by the Dirac distribution $\delta(\ell)$. In turn, the second term accounts for other trajectories that arrived onto $\partial\Omega$ and thus yielded positive $\ell'_\tau$. Note that all the “ingredients” of Eq. (4) depend on the shape of the boundary layer $\Omega_\epsilon$, so that advanced numerical tools are in general needed to access these quantities. In practice, their computation is usually more challenging than that of the boundary local time $\ell_t$ in the original domain $\Omega$. In other words, a direct application of this spectral expansion for Monte Carlo simulations is pointless. In the next subsection, we propose an approximate scheme that relies on the smallness of the width $\epsilon$.

## 2.3. Flat layer approximation (FLA)

For clarity, we focus on the planar case ($d = 2$) but the method is immediately applicable in any dimension $d \ge 2$ (see below). We assume that the boundary $\partial\Omega$ is smooth and can thus be considered as flat at a small scale $\epsilon$. In other words, the boundary layer $\Omega_\epsilon$ in the vicinity of $x_k$ can be approximated by an infinite stripe of width $\epsilon$: $\mathbb{R} \times (0, \epsilon)$. Let $x_b \in \partial\Omega$ be the boundary point, which is the closest to $x_k$. For simplicity of notations, we introduce local coordinates with the origin at $x_b$ and with y axis being perpendicular to the boundary (Fig. 1b). In other words, $x_b = (0, 0)$ and $x_k = (0, y_0)$, where $y_0 = |x_k - \partial\Omega| = |x_k - x_b|$ is the distance to the boundary from $x_k$. As lateral and transverse displacements in the stripe are independent, one can focus on diffusion on the interval $(0, \epsilon)$ with reflecting endpoint 0 and absorbing endpoint $\epsilon$, so that the general spectral expansion can be considerably simplified. In fact, one has

$$j(x, \epsilon, \ell, t|\mathbf{x}_0, y_0) = \frac{e^{-(x-x_0)^2/(4Dt)}}{\sqrt{4\pi Dt}} J(\ell, t|\mathbf{y}_0), \quad (6)$$

where the first factor is the Gaussian probability density of the lateral displacement $x_\tau$, while the second factor describes the joint probability density of $\ell'_\tau$ and $\tau$ for the transverse displacement. According to Eq. (4), the latter admits the following representation in the Laplace domain:

$$\tilde{J}(\ell, p|y_0) = \tilde{J}_0(p|y_0)\delta(\ell) + \sum_{k=0}^{\infty} [V_k^{(p)}(y_0)]^* e^{-\mu_k^{(p)}\ell} C_k^{(p)}, \quad (7)$$

where

$$C_k^{(p)} = -(\partial_n V_k^{(p)})|_{y=\epsilon}, \quad (8)$$

$\mu_k^{(p)}$ and $V_k^{(p)}(y)$ are the eigenvalues and eigenfunctions of the Steklov problem on the interval $(0, \epsilon)$, and $\tilde{J}_0(p|\mathbf{x}_0)$ is the Laplace transform of the probability flux onto the absorbing endpoint $\epsilon$ of the interval under the condition of not hitting the endpoint 0. In the case of an interval $(0, \epsilon)$ with one absorbing endpoint $\epsilon$, the above spectral expansion includes only one eigenmode ($k = 0$), and all the ingredients are known explicitly:

$$V_0^{(p)}(y) = \frac{\sinh(\sqrt{p/D}(\epsilon - y))}{\sinh(\sqrt{p/D}\epsilon)}, \quad (9a)$$

$$\mu_0^{(p)} = \sqrt{p/D} \operatorname{ctanh}(\sqrt{p/D}\epsilon), \quad (9b)$$

$$C_0^{(p)} = \frac{\sqrt{p/D}}{\sinh(\sqrt{p/D}\epsilon)}, \quad (9c)$$

$$\tilde{J}_0(p|y_0) = \frac{\sinh(\sqrt{p/D}y_0)}{\sinh(\sqrt{p/D}\epsilon)}. \quad (9d)$$

As a consequence, one has

$$\tilde{J}(\ell, p|y_0) = \frac{\sinh(\sqrt{p/D}y_0)}{\sinh(\sqrt{p/D}\epsilon)} \delta(\ell) + \frac{\sqrt{p/D} \sinh(\sqrt{p/D}(\epsilon - y_0))}{\sinh^2(\sqrt{p/D}\epsilon)} e^{-\ell} \sqrt{p/D} \operatorname{ctanh}(\sqrt{p/D}\epsilon). \quad (10)$$

In order to generate two random variables $\ell'_\tau$ and $\tau$, one can first generate $\ell'_\tau$ and then employ the conditional distribution of $\tau$. The marginal probability density $\rho(\ell|y_0)$ of the random variable $\ell'_\tau$ is obtained by integrating the joint
probability density $J(\ell, t|y_0)$ over the second variable $t$, i.e., by setting $p = 0$ in the Laplace transform:

$$ \rho(\ell|y_0) = \int_0^\infty dt J(\ell, t|y_0) = \tilde{J}(\ell, 0|y_0) = \pi_0 \delta(\ell) + (1 - \pi_0) \frac{e^{-\ell/\epsilon}}{\epsilon}, \quad (11) $$

where $\pi_0 = \tilde{J}_0(0|y_0) = y_0/\epsilon$ is the splitting probability of reaching the endpoint $\epsilon$ first. This is a mixture of the Dirac distribution at 0 and the exponential distribution with mean $\epsilon$. In other words, the random variable $\ell'_\tau$ admits a simple form:

$$ \ell'_{\tau} = \begin{cases} 0, & \text{with probability } \pi_0, \\ \epsilon\eta, & \text{with probability } 1 - \pi_0, \end{cases} \quad (12) $$

where $\eta$ is the standard exponentially distributed random variable (with mean 1).

Even though the conditional distribution of the escape time $\tau$ is found exactly in the Laplace domain (see below), the inversion of the Laplace transform and the consequent generation of $\tau$ can be time-consuming (see Appendix B). To overcome this limitation, we propose a simple approximation, in which the random escape time $\tau$ is replaced by its mean value, which can be found explicitly (see Appendix C for discussion on its validity).

In practice, one can first generate a binary random variable that determines which endpoint of the interval $(0, \epsilon)$ is reached first: the endpoint $\epsilon$ (with probability $\pi_0$) or the endpoint 0 (with probability $1 - \pi_0$).

(i) If the endpoint $\epsilon$ is reached first, then the boundary local time $\ell'_\tau$ is zero, while the escape time (denoted here as $\tau_0$) is determined by the conditional probability density $H_0(t|y_0)$ (see Appendix B), whose Laplace transform is

$$ \langle e^{-p\tau_0} \rangle = \tilde{H}_0(p|y_0) = \frac{\epsilon}{y_0} \frac{\sinh(\sqrt{p/D} y_0)}{\sinh(\sqrt{p/D} \epsilon)}. \quad (13) $$

The mean escape time is then

$$ \langle \tau_0 \rangle = - \left( \partial_p \tilde{H}_0(p|y_0) \right)_{p=0} = \frac{\epsilon^2 - y_0^2}{6D}. \quad (14) $$

(ii) If the endpoint 0 is reached first, then the boundary local time $\ell'_\tau$ is nonzero and can be generated from its marginal distribution (11), i.e., by setting $\ell'_\tau = \epsilon\eta$, where $\eta$ is the standard exponential random variable. In turn, the Laplace transform of the conditional probability density of the escape time (denoted here as $\tau_\ell$) is determined by the second term in Eq. (10), which is divided by the second term in Eq. (11):

$$ \langle e^{-p\tau_\ell} \rangle = \tilde{H}(p|y_0, \ell) = \frac{1}{(1 - \pi_0)e^{-\ell/\epsilon}/\epsilon} \frac{\sqrt{p/D} \sinh(\sqrt{p/D}(\epsilon - y_0))}{\sinh^2(\sqrt{p/D}\epsilon)} e^{-\ell} \sqrt{p/D} \tanh(\sqrt{p/D}\epsilon). \quad (15) $$

The mean escape time is then

$$ \langle \tau_\ell \rangle = - \left( \partial_p \tilde{H}(p|y_0, \ell) \right)_{p=0} = \frac{2\ell\epsilon + 2\epsilon^2 - (\epsilon - y_0)^2}{6D}. \quad (16) $$

Finally, the lateral displacement $x_\tau$ can be generated as the Gaussian random variable with mean zero and variance $2D\langle\tau_0\rangle$ or $2D\langle\tau_\ell\rangle$.

In higher dimensions ($d > 2$), the boundary layer near $x_k$ is approximated by a slab $\mathbb{R}^{d-1} \times (0, \epsilon)$ between two parallel (hyper)planes. As the lateral and transverse displacements are still independent, the probability density flux $j(x, \ell, t|x_k)$ is again factored into a Gaussian density for the lateral displacement and the same joint probability density $J(\ell, t|y_0)$ for the interval $(0, \epsilon)$, like in Eq. (6). As a consequence, the generation of $\ell'_\tau$ and $\tau$ remains unchanged, whereas the displacement along each lateral coordinate is still an independent Gaussian random variable with mean 0 and variance $2D\langle\tau_0\rangle$ or $2D\langle\tau_\ell\rangle$.

## 2.4. Curved layer approximation (CLA)

The above flat layer approximation ignores the curvature of the boundary, that is only valid for small enough $\epsilon$. Such a choice may slow down simulations, as the process stays too long in the vicinity of the boundary. This approximation can be further improved by approximating the boundary layer by a circular annulus or a spherical shell, with the same radius of curvature. We start by the two-dimensional case and then extend it to three dimensions.
<figure><img src="image_3.png" /><figcaption>Figure 2: Curved layer approximation: starting from a point $(r_0, \phi_0)$ in polar coordinates (yellow circle), a random trajectory (in blue) experiences numerous reflections on the circle of radius $R$ before escaping a thin layer of width $\epsilon$ at a random escape time $\tau$ in a random escape position $(R \pm \epsilon, \phi_\tau)$ (red circle), with sign minus corresponding to a convex boundary (a) and sign plus corresponding to a concave boundary (b). These reflections are characterized by the acquired boundary local time $\ell'_\tau$.</figcaption></figure>

### 2.4.1. Convex boundary in two dimensions

We first consider the case of a convex boundary (Fig. 2a). Let $R$ be the radius of curvature of the boundary at the boundary point $x_b$, which is the closest to $x_k$. The boundary layer can be approximated by a circular annulus between two circles of radii $R - \epsilon$ and $R$, separated by distance $\epsilon$, i.e., $\Omega_\epsilon = \{x \in \mathbb{R}^2 : R - \epsilon < |x| < R\}$. We are interested in the escape event from that annulus through the inner circle $\Gamma_\epsilon = \{x \in \mathbb{R}^2 : |x| = R - \epsilon\}$ of radius $R - \epsilon$, which is thus absorbing. The joint probability density $J(\ell, t|x_k)$ of the escape time $\tau$ and the boundary local time $\ell'_\tau$, acquired on the outer circle up to $\tau$, is determined by the radial displacements. This can also be seen by integrating the probability density flux $j(x, \ell, t|x_k)$ over the escape position $X'_\tau$ on the inner circle, in which case all terms in the sum of Eq. (4) with $k \neq 0$ vanish due to the periodic dependence of the Steklov eigenfunctions $V_k^{(p)}$ on the angle [24]. In other words, we get

$$ \tilde{J}(\ell, p|r_0) = \tilde{J}_0(p|r_0)\delta(\ell) + V_0^{(p)}(r_0)C_0^{(p)}e^{-\mu_0^{(p)}\ell}, \quad (17) $$

where $r_0 = |x_k|$ is the radial coordinate of the starting point (here, we do not employ local coordinates centered at $x_b$ so that $r_0 \in (R - \epsilon, R])$, and

$$ C_0^{(p)} = - \int_{\Gamma_\epsilon} d\mathbf{x} (\partial_n V_0^{(p)}). \quad (18) $$

For a circular annulus with an absorbing inner boundary, one has [24]

$$ V_0^{(p)}(r_0) = \frac{1}{\sqrt{2\pi R}} \frac{K_0(\alpha(R - \epsilon))I_0(\alpha r_0) - I_0(\alpha(R - \epsilon))K_0(\alpha r_0)}{K_0(\alpha(R - \epsilon))I_0(\alpha R) - I_0(\alpha(R - \epsilon))K_0(\alpha R)}, \quad (19a) $$

$$ \mu_0^{(p)} = \alpha \frac{K_0(\alpha(R - \epsilon))I_1(\alpha R) + I_0(\alpha(R - \epsilon))K_1(\alpha R)}{K_0(\alpha(R - \epsilon))I_0(\alpha R) - I_0(\alpha(R - \epsilon))K_0(\alpha R)}, \quad (19b) $$

$$ C_0^{(p)} = \frac{\sqrt{2\pi}}{\sqrt{R}} \frac{1}{K_0(\alpha(R - \epsilon))I_0(\alpha R) - I_0(\alpha(R - \epsilon))K_0(\alpha R)}, \quad (19c) $$

$$ \tilde{J}_0(p|r_0) = \frac{K_0(\alpha R)I_0(\alpha r_0) - I_0(\alpha R)K_0(\alpha r_0)}{K_0(\alpha R)I_0(\alpha(R - \epsilon)) - I_0(\alpha R)K_0(\alpha(R - \epsilon))}, \quad (19d) $$

where $\alpha = \sqrt{p/D}$, $I_\nu(z)$ and $K_\nu(z)$ are the modified Bessel functions of the first and second kind, and we used $I'_0(z) = I_1(z)$, $K'_0(z) = -K_1(z)$, and the Wronskian $K_\nu(z)I'_\nu(z) - I_\nu(z)K'_\nu(z) = 1/z$.

In the limit $p \to 0$, we get the marginal probability density of the boundary local time $\ell'_\tau$:

$$ \rho(\ell|r_0) = \tilde{J}(\ell, 0|r_0) = \tilde{J}_0(0|r_0)\delta(\ell) + V_0^{(0)}(r_0)C_0^{(0)}e^{-\mu_0^{(0)}\ell}, \quad (20) $$
where

$$V_0^{(0)}(r_0) = \frac{1}{\sqrt{2\pi R}} \frac{\ln(r_0/(R-\epsilon))}{\ln(R/(R-\epsilon))}, \qquad (21a)$$

$$\mu_0^{(0)} = \frac{1}{R \ln(R/(R-\epsilon))}, \qquad (21b)$$

$$C_0^{(0)} = \frac{\sqrt{2\pi}}{\sqrt{R}} \frac{1}{\ln(R/(R-\epsilon))}, \qquad (21c)$$

$$\pi_0 = \tilde{J}_0(0|r_0) = \frac{\ln(R/r_0)}{\ln(R/(R-\epsilon))}, \qquad (21d)$$

and we used the asymptotic behavior of the modified Bessel functions. As a consequence, we find

$$\rho(\ell|r_0) = \pi_0 \delta(\ell) + (1 - \pi_0) \mu_0^{(0)} e^{-\mu_0^{(0)} \ell}. \qquad (22)$$

As in the case of the flat layer approximation, we get a mixture of the Dirac distribution at 0 and the exponential distribution with the mean $1/\mu_0^{(0)}$, where $\pi_0$ is the splitting probability of the first arrival onto the absorbing (inner) circle, in which case the boundary local time is zero. Setting $r_0 = R - y_0$, one sees that $\pi_0(r_0) \approx y_0/\epsilon$ in the limit $\epsilon \ll 1$ and $y_0 \ll 1$, as in the flat case.

By generating a binary random variable, one can choose between two options:

(i) with probability $\pi_0$, one has $\ell'_\tau = 0$, and the escape time $\tau_0$ is determined by

$$\langle e^{-p\tau_0} \rangle = \tilde{H}_0(p|r_0) = \frac{\tilde{J}_0(p|r_0)}{\pi_0}. \qquad (23)$$

In particular, the mean escape time is

$$\langle \tau_0 \rangle = - \left( \partial_p \tilde{H}_0(p|r_0) \right)_{p=0} = \frac{1}{4D} \left( (R-\epsilon)^2 - r_0^2 + \frac{R^2 - r_0^2}{\ln(R/r_0)} - \frac{R^2 - (R-\epsilon)^2}{\ln(R/(R-\epsilon))} \right). \qquad (24)$$

(ii) with probability $1 - \pi_0$, one has $\ell'_\tau = \eta/\mu_0^{(0)}$ (where $\eta$ is the standard exponential random variable), while the escape time $\tau_\ell$ is determined by

$$\langle e^{-p\tau_\ell} \rangle = \tilde{H}(p|r_0, \ell) = \frac{V_0^{(p)}(r_0) C_0^{(p)} e^{-\mu_0^{(p)} \ell}}{(1 - \pi_0) \mu_0^{(0)} e^{-\mu_0^{(0)} \ell}}. \qquad (25)$$

The mean escape time is then

$$\langle \tau_\ell \rangle = - \left( \partial_p \tilde{H}(p|r_0, \ell) \right)_{p=0} = \frac{\ell V_0^{(0)} C_0^{(0)} (\partial_p \mu_0^{(p)})_{p=0} - (\partial_p V_0^{(p)}(r_0) C_0^{(p)})_{p=0}}{(1 - \pi_0) \mu_0^{(0)}}, \qquad (26)$$

where

$$(\partial_p \mu_0^{(p)})_{p=0} = R \frac{1 - (1 - \varepsilon/R)^2 + 2 \ln(1 - \varepsilon/R) (\ln(1 - \varepsilon/R) + 1)}{4D \ln^2(1 - \varepsilon/R)}, \qquad (27)$$

and

$$(\partial_p V_0^{(p)}(r_0) C_0^{(p)})_{p=0} = - \frac{1}{4DR \ln^3(R/(R-\epsilon))} \left( (2R^2 + (R-\epsilon)^2 - r_0^2) \ln(r_0/(R-\epsilon)) \ln(R/(R-\epsilon)) + 2(R^2 - (R-\epsilon)^2) \ln(R/r_0) - (2R^2 - (R-\epsilon)^2 - r_0^2) \ln(R/(R-\epsilon)) \right). \qquad (28)$$

Substituting these expressions into Eq. (26), one can evaluate the mean escape time $\langle \tau_\ell \rangle$.

The last step consists in generating the escape position $X'_\tau$ on the inner circle. In contrast to the flat case, “lateral” displacements along the tangential direction are not independent from the transverse (radial) displacements due to the
curvature. The exact distribution of the escape position can be found from the probability density flux $j(\mathbf{x}, \ell, t|\mathbf{x}_0)$, but its practical computation could be too sophisticated. In turn, two approximations can be useful.

(i) The simplest approximation consists in replacing the random variable $X'_\tau$ by its mean value. By symmetry, one has $\langle X'_\tau \rangle = (R - \varepsilon, \phi_0)$ in polar coordinates, where $\phi_0$ is the angle of the starting position $\mathbf{x}_k = (r_0, \phi_0)$. For a flat boundary, this approximation would correspond to replacing $x_\tau$ by $\langle x_\tau \rangle = 0$.

(ii) This approximation can be improved for small $\varepsilon$, for which the lateral motion can be treated as lying on the inner circle of radius $R - \varepsilon$ (instead of displacements inside the thin annulus of width $\varepsilon$). In this approximation, the displacements on the circle are characterized by the angle $\phi_\ell$, which is independent of radial displacements, and can thus be generated from known distributions (see [61, 62]).

### 2.4.2. Concave boundary in two dimensions

The analysis for a locally concave boundary (Fig. 2b) is very similar. In this case, one considers a thin circular annulus $\Omega_\varepsilon = \{\mathbf{x} \in \mathbb{R}^2 : R < |\mathbf{x}| < R + \varepsilon\}$ with the absorbing outer boundary $\Gamma_\varepsilon = \{\mathbf{x} \in \mathbb{R}^2 : |\mathbf{x}| = R + \varepsilon\}$. As the absorbing boundary is now located at $R + \varepsilon$ instead of $R - \varepsilon$, one mainly needs to change the sign of $\varepsilon$:

$$V_0^{(p)}(r_0) = \frac{1}{\sqrt{2\pi R}} \frac{K_0(\alpha(R+\varepsilon))I_0(\alpha r_0) - I_0(\alpha(R+\varepsilon))K_0(\alpha r_0)}{K_0(\alpha(R+\varepsilon))I_0(\alpha R) - I_0(\alpha(R+\varepsilon))K_0(\alpha R)}, \quad (29a)$$

$$\mu_0^{(p)} = -\alpha \frac{K_0(\alpha(R+\varepsilon))I_1(\alpha R) + I_0(\alpha(R+\varepsilon))K_1(\alpha R)}{K_0(\alpha(R+\varepsilon))I_0(\alpha R) - I_0(\alpha(R+\varepsilon))K_0(\alpha R)}, \quad (29b)$$

$$C_0^{(p)} = -\frac{\sqrt{2\pi}}{\sqrt{R}} \frac{1}{K_0(\alpha(R+\varepsilon))I_0(\alpha R) - I_0(\alpha(R+\varepsilon))K_0(\alpha R)}, \quad (29c)$$

$$\tilde{J}_0(p|r_0) = \frac{K_0(\alpha R)I_0(\alpha r_0) - I_0(\alpha R)K_0(\alpha r_0)}{K_0(\alpha R)I_0(\alpha(R+\varepsilon)) - I_0(\alpha R)K_0(\alpha(R+\varepsilon))}. \quad (29d)$$

Note also that the signs in front of $\mu_0^{(p)}$ and $C_0^{(p)}$ are also changed due to the opposite direction of the normal derivative. These expressions determine $\tilde{J}(\ell, p|r_0)$ via Eq. (17). One also gets

$$V_0^{(0)}(r_0) = \frac{1}{\sqrt{2\pi R}} \frac{\ln(r_0/(R+\varepsilon))}{\ln(R/(R+\varepsilon))}, \quad (30a)$$

$$\mu_0^{(0)} = -\frac{1}{R \ln(R/(R+\varepsilon))}, \quad (30b)$$

$$C_0^{(0)} = -\frac{\sqrt{2\pi}}{\sqrt{R}} \frac{1}{\ln(R/(R+\varepsilon))}, \quad (30c)$$

$$\pi_0 = \tilde{J}_0(0|r_0) = \frac{\ln(R/r_0)}{\ln(R/(R+\varepsilon))}, \quad (30d)$$

that determine $\rho(\ell|r_0)$ via Eq. (20). Finally, the mean escape time is given by

$$\langle \tau_0 \rangle = \frac{1}{4D} \left( (R + \varepsilon)^2 - r_0^2 + \frac{R^2 - r_0^2}{\ln(R/r_0)} - \frac{(R + \varepsilon)^2 - R^2}{\ln(1 + \varepsilon/R)} \right), \quad (31)$$

whereas $\langle \tau_\ell \rangle$ is still determined by Eq. (26), in which

$$(\partial_p \mu_0^{(p)})_{p=0} = -R \frac{1 - (1 + \varepsilon/R)^2 + 2 \ln(1 + \varepsilon/R)(\ln(1 + \varepsilon/R) + 1)}{4D \ln^2(1 + \varepsilon/R)}, \quad (32)$$

and

$$(\partial_p V_0^{(p)}(r_0)C_0^{(p)})_{p=0} = \frac{1}{4DR \ln^3(R/(R+\varepsilon))} \left( (2R^2 + (R+\varepsilon)^2 - r_0^2) \ln(r_0/(R+\varepsilon)) \ln(R/(R+\varepsilon)) + 2(R^2 - (R+\varepsilon)^2) \ln(R/r_0) - (2R^2 - (R+\varepsilon)^2 - r_0^2) \ln(R/(R+\varepsilon)) \right). \quad (33)$$
### 2.4.3. Convex boundary in three dimensions

A natural extension of the above analysis to the three-dimensional case would consist in approximating a smooth boundary with two radii of curvature by an ellipsoidal or hyperboloidal layer. Unfortunately, there is no explicit representation of the joint probability density of $l'_τ$ and $\tau$ in these cases. We therefore restrict the discussion to the particular setting, in which the smooth boundary can be locally approximated by a sphere of radius $R$ (i.e. two radii of curvature are equal). In this case, the extension is straightforward, while the resulting formulas are actually simpler than in two dimensions.

A convex boundary is approximated by a spherical shell $\Omega_\varepsilon = \{x \in \mathbb{R}^3 : R - \varepsilon < |x| < R\}$ with the absorbing inner sphere $\Gamma_\varepsilon = \{x \in \mathbb{R}^3 : |x| = R - \varepsilon\}$. In this case, the modified Bessel functions $I_\nu(z)$ and $K_\nu(z)$ are replaced by the modified spherical Bessel functions $i_\nu(z)$ and $k_\nu(z)$, which admit simple expressions. Skipping technical steps, we get

$$V_0^{(p)}(r_0) = \frac{\sinh(\alpha(r_0 - R + \varepsilon))}{\sqrt{4\pi} r_0 \sinh(\alpha\varepsilon)}, \qquad (34a)$$

$$\mu_0^{(p)} = \alpha \tanh(\alpha\varepsilon) - \frac{1}{R}, \qquad (34b)$$

$$C_0^{(p)} = \frac{\sqrt{4\pi}\alpha(R - \varepsilon)}{\sinh(\alpha\varepsilon)}, \qquad (34c)$$

$$\tilde{J}_0(p|r_0) = \frac{(R - \varepsilon) \sinh(\alpha(R - r_0))}{r_0 \sinh(\alpha\varepsilon)}, \qquad (34d)$$

while their limits as $p \to 0$ are

$$V_0^{(0)}(r_0) = \frac{r_0 - R + \varepsilon}{\sqrt{4\pi}r_0\varepsilon}, \qquad (35a)$$

$$\mu_0^{(0)} = \frac{1}{\varepsilon} - \frac{1}{R}, \qquad (35b)$$

$$C_0^{(0)} = \frac{\sqrt{4\pi(R - \varepsilon)}}{\varepsilon}, \qquad (35c)$$

$$\pi_0 = \tilde{J}_0(0|r_0) = \frac{(R - \varepsilon)(R - r_0)}{\varepsilon r_0}. \qquad (35d)$$

As a consequence, one can easily evaluate the mean escape times:

$$\langle \tau_0 \rangle = \frac{\varepsilon^2 - (R - r_0)^2}{6D}, \qquad (36a)$$

$$\langle \tau_\ell \rangle = \frac{2\ell\varepsilon + 2\varepsilon^2 - (r_0 - R + \varepsilon)^2}{6D}. \qquad (36b)$$

Setting $y_0 = R - r_0$ as the distance to the boundary, one retrieves Eqs. (14, 16) for the flat case. As previously, one first employs the splitting probability $\pi_0$ to choose between two options, generates the boundary local time $l'_\tau$, and uses the mean escape time. The escape position can be either generated, or replaced by its mean value, which by symmetry is $(R - \varepsilon, \theta_0, \phi_0)$ in spherical coordinates, where $\theta_0$ and $\phi_0$ are the angular coordinates of the starting position $x_k$.
### 2.4.4. Concave boundary in three dimensions

A concave boundary is approximated by a spherical shell $\Omega_{\varepsilon} = \{\mathbf{x} \in \mathbb{R}^3 : R < |\mathbf{x}| < R + \varepsilon\}$ with the absorbing outer sphere $\Gamma_{\varepsilon} = \{\mathbf{x} \in \mathbb{R}^3 : |\mathbf{x}| = R + \varepsilon\}$. One has

$$V_0^{(p)}(r_0) = \frac{\sinh(\alpha(R - r_0 + \varepsilon))}{\sqrt{4\pi} r_0 \sinh(\alpha\varepsilon)}, \qquad (37a)$$

$$\mu_0^{(p)} = \alpha \tanh(\alpha\varepsilon) + \frac{1}{R}, \qquad (37b)$$

$$C_0^{(p)} = \frac{\sqrt{4\pi}\alpha(R + \varepsilon)}{\sinh(\alpha\varepsilon)}, \qquad (37c)$$

$$\tilde{J}_0(p|r_0) = \frac{(R + \varepsilon) \sinh(\alpha(r_0 - R))}{r_0 \sinh(\alpha\varepsilon)}, \qquad (37d)$$

and

$$V_0^{(0)}(r_0) = \frac{R - r_0 + \varepsilon}{\sqrt{4\pi r_0 \varepsilon}}, \qquad (38a)$$

$$\mu_0^{(0)} = \frac{1}{\varepsilon} + \frac{1}{R}, \qquad (38b)$$

$$C_0^{(0)} = \frac{\sqrt{4\pi(R + \varepsilon)}}{\varepsilon}, \qquad (38c)$$

$$\pi_0 = \tilde{J}_0(0|r_0) = \frac{(R + \varepsilon)(r_0 - R)}{\varepsilon r_0}. \qquad (38d)$$

As a consequence, one can easily evaluate the mean escape times:

$$\langle \tau_0 \rangle = \frac{\varepsilon^2 - (r_0 - R)^2}{6D}, \qquad (39a)$$

$$\langle \tau_{\ell} \rangle = \frac{2\ell\varepsilon + 2\varepsilon^2 - (R - r_0 + \varepsilon)^2}{6D}. \qquad (39b)$$

Setting $y_0 = r_0 - R$ as the distance to the boundary, one retrieves Eqs. (14, 16) for the flat case.

### 2.5. Stopping conditions

In the conventional setting of domains with an absorbing boundary, the simulation was stopped upon the first arrival into the boundary layer of width $\varepsilon$. In turn, as we deal with reflected Brownian motion, the simulation would never stop without imposing a suitable stopping condition. In many applications, one employs one of the three conditions (or their combinations):

(i) at a fixed time $T$;

(ii) at a random stopping time $\delta$ with a prescribed probability distribution;

(iii) at a random time $\mathcal{T}$ when the boundary local time $\ell_t$ exceeds a fixed or random threshold $\hat{\ell}$ with a prescribed probability distribution.

The first stopping condition allows one to characterize the diffusive dynamics in a confining domain $\Omega$, i.e., the position $X_T$ and the boundary local time $\ell_T$ at a given time $T$. In particular, one can obtain the empirical distributions of these quantities, their moments and correlations. The second stopping condition is often employed to simulate “mortal” particles with a finite random lifetime with a given distribution [63–65]. The most common choice is the exponential distribution of the stopping time $\delta$, $\mathbb{P}\{\delta > t\} = e^{-pt}$, with $p$ being the decay rate or, equivalently, $1/p$ being the mean lifetime. Basic examples are nuclear radioactive disintegration, spin relaxation, photobleaching, first-order bulk reactions, etc. In this case, one generates the statistics of the position $X_\delta$ and the acquired boundary local time $\ell_\delta$ until the “death” time $\delta$. The third stopping condition was introduced in [4, 52] to describe various surface reactions which may occur after a sufficient number (represented by the threshold $\hat{\ell}$) of failed reaction attempts. In this way, one
can generate the statistics of the reaction times $\mathcal{T}$ and of the reaction locations $X_T$ on the boundary. If the threshold $\hat{l}$ is fixed, the surface reaction is triggered after $\hat{l}/\epsilon$ encounters of the diffusing particle with the reactive layer of width $\epsilon$. More generally, the number of required encounters can be drawn from a given distribution. For instance, the most common setting of a constant partial reactivity $\kappa$ in the Robin boundary condition corresponds to the exponential distribution of the threshold, $\mathbb{P}\{\hat{l} > \ell\} = e^{-q\ell}$, with $q = \kappa/D$ [4]. In turn, other distributions can account for more sophisticated surface reactions such as activation or passivation of the reactive boundary, or non-Markovian binding [15, 16]. The third stopping condition was also applied to describe permeation across the boundary [17–19].

An accurate implementation of the stopping condition is known to be delicate even for the absorbing boundary, i.e., for the classical WOS algorithm. For instance, if the stopping time $T$ is fixed, the simulation does not actually stop at $T$ but when $t_{k+1} > T$. However, if the last jump between $t_k$ and $t_{k+1}$ was big, the position $X_T$ lying inside the large ball $B_{\rho_k}(x_k)$ can significantly differ from the last generated position $X_{t_{k+1}}$ on the sphere $\partial B_{\rho_k}(x_k)$. To amend this problem, one often employs specific treatment of the last step such as, e.g., generating the random position $X_T$ for the intermediate time $T$ between $t_k$ and $t_{k+1}$ from the conditional probability density, an interpolation between $X_{t_k}$ and $X_{t_{k+1}}$, the inclusion of the upper bound on the possible escape time, etc [66, 67]. As this issue is classical, we do not discuss these improvements. In our work, we focus on the statistics of the boundary local time, for which the last step is less critical. Indeed, if the last step (after which the simulation is stopped) is a big jump to the sphere $\partial B_{\rho_k}(x_k)$, the bulk trajectory inside the ball $B_{\rho_k}(x_k)$ cannot touch the boundary, so that the boundary local time has not changed during this jump: $\ell_{t_{k+1}} = \ell_T = \ell_{t_k}$. In turn, if the last jump is an escape from the boundary layer $\Omega_\epsilon$ of width $\epsilon$, the boundary local time can change, but its random increment $\ell'_\tau$ is typically of the order of $\epsilon$. If the width $\epsilon$ is small enough, the induced error of $\ell_T$ is small due to $\ell_{t_k} \le \ell_T \le \ell_{t_{k+1}}$. As a consequence, either of these two bounds can be used to approximate $\ell_T$ (moreover, one can record both $\ell_{t_k}$ and $\ell_{t_{k+1}}$ to get lower and upper bounds). Alternatively, one can employ a suitable interpolation between $\ell_{t_k}$ and $\ell_{t_{k+1}}$, e.g., the Skorokhod integral representation of the boundary local time (see Appendix D) suggests to set

$$ \ell_T = \ell_{t_k} + (\ell_{t_{k+1}} - \ell_{t_k}) \frac{\sqrt{T} - \sqrt{t_k}}{\sqrt{t_{k+1}} - \sqrt{t_k}}. \quad (40) $$

Finally, one can also attempt to derive the conditional probability density of $\ell_T$. In the following numerical examples, we employ the simplest approximation by setting $\ell_T = \ell_{t_{k+1}}$.

### 3. Implementation of the algorithm

We summarize the main steps of our method by focusing explicitly on domains in two and three dimensions. An extension to higher dimensions is immediately available for FLA; in turn, additional computations would be needed for implementing CLA (not present).

1. Consider a domain $\Omega$ defined by a function $\rho(x)$, which determines the distance from any point $x$ to the boundary $\partial\Omega$. Initialize all the parameters: initial position $x_{\text{ini}}$, diffusion coefficient $D$, number of particles $N$, width $\epsilon$ of the boundary layer, stopping condition, etc. Load the pre-computed data to generate the escape times for the unit disk in 2d, or the unit sphere in 3d (see Appendix A). For each particle ($i = 1, 2, \dots, N$), run the following simulation.

2. Initiate $x_0 = x_{\text{ini}}$, $t_0 = 0$, $\ell_0 = 0$; for the second stopping condition, generate the random stopping time $\delta$ (e.g., if $\delta$ obeys the exponentially distribution, set $\delta = -\ln(u)/\rho$, where $u$ is the uniform random variable on (0, 1), and $p > 0$ is a given parameter); for the third stopping condition, generate the random threshold $\hat{l}$ (e.g., if $\hat{l}$ obeys the exponentially distribution, set $\hat{l} = -\ln(u)/q$, where $u$ is the uniform random variable on (0, 1), and $q > 0$ is a given parameter).

3. At each step $k = 0, 1, 2, \dots$, the radius $\rho_k = \rho(x_k)$ is determined by the current position $x_k$.

* If $\rho_k > \epsilon/2$, the new position $x_{k+1}$ is generated to be uniformly distributed on $\partial B_{\rho_k}(x_k)$. In two dimensions, it is the circle of radius $\rho_k$, centered at $x_k$, so that

$$ \begin{cases} x_{k+1} = x_k + \rho_k \cos \varphi, \\ y_{k+1} = y_k + \rho_k \sin \varphi, \end{cases} \quad (41) $$
<table><caption>Table 1: Summary of approximations for the quantities that determine the escape-from-a-layer event.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">FLA</th><th colspan="4">CLA</th></tr><tr><th>convex 2d</th><th>concave 2d</th><th>convex 3d</th><th>concave 3d</th></tr></thead><tbody><tr><th>π<sub>0</sub></th><td>Eq. (11)</td><td>Eq. (21d)</td><td>Eq. (30d)</td><td>Eq. (35d)</td><td>Eq. (38d)</td></tr><tr><th>⟨τ<sub>0</sub>⟩</th><td>Eq. (14)</td><td>Eq. (24)</td><td>Eq. (31)</td><td>Eq. (36a)</td><td>Eq. (39a)</td></tr><tr><th>⟨τ<sub>ℓ</sub>⟩</th><td>Eq. (16)</td><td>Eq. (26)</td><td>Eq. (26)</td><td>Eq. (36b)</td><td>Eq. (39b)</td></tr><tr><th>μ<sub>0</sub><sup>(0)</sup></th><td>Eq. (12)</td><td>Eq. (21b)</td><td>Eq. (30b)</td><td>Eq. (35b)</td><td>Eq. (38b)</td></tr></tbody></table>

where $\varphi$ is a random variable uniformly distributed in $(0, 2\pi)$. In three dimensions, it is the sphere of radius $\rho_k$, centered at $\mathbf{x}_k$, so that

$$ \begin{cases} x_{k+1} = x_k + \rho_k \sin \theta \cos \varphi, \\ y_{k+1} = y_k + \rho_k \sin \theta \sin \varphi, \\ z_{k+1} = z_k + \rho_k \cos \theta, \end{cases} \qquad (42) $$

where $\varphi$ is a random variable uniformly distributed in $(0, 2\pi)$, and $\theta$ is a random variable such that $\cos \theta$ is uniformly distributed in $(-1, 1)$.

The diffusion time is incremented by

$$ t_{k+1} = t_k + \frac{\rho_k^2}{D} \tau_c, \qquad (43) $$

where $\tau_c$ is the escape time either from the unit disk (in 2d) or from the unit ball (in 3d). This random variable is generated in a standard way (see Appendix A).

* If $\rho < \varepsilon/2$, initiate the escape-from-a-layer event: generate a random variable $w$ uniformly distributed in $(0, 1)$.

- If $w \in (0, \pi_0)$, the particle escapes the layer without touching the boundary, so that $\ell'_{\tau} = 0$. We set therefore $t_{k+1} = t_k + \langle \tau_0 \rangle$ and $\ell_{t_{k+1}} = \ell_{t_k}$.

- If $w \in (\pi_0, 1)$, the particle escapes the layer after reflections on the boundary. The boundary local time increases as $\ell_{t_{k+1}} = \ell_{t_k} + \eta/\mu_0^{(0)}$, where $\eta$ is a standard exponentially random variable, and $t_{k+1} = t_k + \langle \tau_\ell \rangle$.

According to different approximations (FLA or CLA), the curvature sign (convex or concave), and space dimension ($d=2$ or $d=3$), $\pi_0$, $\langle \tau_0 \rangle$, $\langle \tau_\ell \rangle$, and $\mu_0^{(0)}$ are given by equations listed in Table 1, with $R$ being the radius of curvature at the boundary point $\mathbf{x}_b$ closest to $\mathbf{x}_k$, and $r_0 = |R - \rho_k|$. The new position $\mathbf{x}_{k+1}$ is located at distance $\varepsilon$ from the boundary point $\mathbf{x}_b$.

4. Check for the stopping condition ($t_{k+1} > T$, or $t_{k+1} > \delta$, or $\ell_{t_{k+1}} > \hat{\ell}$): if it is not satisfied, return to step 3; otherwise, stop the simulation and record the boundary local time $\ell_{t_{k+1}}$ for this particle.

5. Repeat the steps 2-4 for $N$ particles to get an empirical statistics of boundary local times.

In the above implementation, the escape-from-a-layer event is triggered when the distance to the boundary $\rho_k$ is below $\varepsilon/2$ (instead of the earlier discussed condition $\rho_k < \varepsilon$). This “trick” allows one to eliminate too often immediate returns to the layer after the escape. The choice of the fraction $1/2$ in front of $\varepsilon$ is somewhat arbitrary, it can be replaced by any number between 0 and 1.

## 4. Validation

In order to check the accuracy of the proposed algorithm and the role of its parameters (such as $\varepsilon$), we consider reflected Brownian motion in a bounded domain $\Omega \subset \mathbb{R}^d$ started from a fixed point $\mathbf{x}_0 \in \Omega$. The probability density $\rho(\ell, t|\mathbf{x}_0)$ of the boundary local time $\ell_t$ reads in the Laplace domain as [23]

$$ \tilde{\rho}(\ell, p|\mathbf{x}_0) = \int_0^\infty dt e^{-pt} \rho(\ell, t|\mathbf{x}_0) = \tilde{S}_0(p|\mathbf{x}_0)\delta(\ell) + \sum_{k=0}^\infty e^{-\mu_k^{(p)}\ell} c_k^{(p)}(\mathbf{x}_0), \qquad (44) $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Considered examples. (a) Unit disk centered at the origin; the initial point (in red) is located at ($r_0 = 0.5, \varphi_0 = 0$). (b) Annulus between two concentric circles of radii $R = 1, L = 2$; the initial point (in red) is located at ($r_0 = 1.5, \varphi_0 = 0$). (c) Unit disk ($R = 1$) centered at the origin and perforated by ten small circular holes of the same radius $R_0 = 0.1$ located randomly (their centers are at $(-0.4167, 0.4213), (0.0973, 0.6397), (0.6418, 0.4860), (-0.7173, 0.0827), (-0.2075, 0.0666), (0.1105, -0.0605), (0.6176, 0.1071), (-0.5381, -0.4460), (-0.0731, -0.5457), (0.3528, -0.5425)$); the initial points are: ($R, 0), (0, 0$). (d) Unit sphere ($R = 1$) centered at the origin; the initial point (in red) is located at $r_0 = 0.5$. (e) Spherical shell between two concentric spheres of radii $R = 1$ and $L = 2$; the initial point (in red) is located at $r_0 = 1.5$. (f) Torus with $R = 1.5, R_0 = 0.5$; the initial point (in red) is located at $(1.5, 0, 0)$.</figcaption></figure>

where

$$c_k^{(p)}(\mathbf{x}_0) = \frac{[V_k^{(p)}(\mathbf{x}_0)]^*}{D} \int_{\Omega} d\mathbf{x} V_k^{(p)}(\mathbf{x}), \quad (45)$$

and $\tilde{S}_0(p|\mathbf{x}_0)$ is the Laplace-transformed survival probability in the presence of an absorbing boundary $\partial\Omega$. While $\tilde{S}_0(p|\mathbf{x}_0)$, $\mu_k^{(p)}$ and $V_k^{(p)}(\mathbf{x})$ can be found numerically for a given domain $\Omega$ (e.g., by a finite-element method), the need for the Laplace transform inversion makes the above spectral expansion less suitable for validating Monte Carlo simulations.

To overcome this difficulty, we focus on the boundary local time $\ell_\delta$, stopped at an exponentially distributed random time $\delta$, i.e., $\mathbb{P}\{\delta > t\} = e^{-pt}$, with a given rate $p > 0$ (the second stopping condition). The probability density of $\ell_\delta$ can be obtained by averaging $\ell_t$ as

$$\Psi(\ell, p|\mathbf{x}_0) = \int_0^\infty dt \underbrace{p e^{-pt}}_{= \text{pdf of } \delta} \rho(\ell, t|\mathbf{x}_0) = p \tilde{\rho}(\ell, p|\mathbf{x}_0). \quad (46)$$

In other words, the statistics of $\ell_\delta$ is directly accessible in the Laplace domain and can thus be obtained without Laplace transform inversion.

We test our approach for six domains shown in Fig. 3. For these domains, we obtain the empirical statistics of $\ell_\delta$ by using FLA and CLA with different $\varepsilon$ and compare to either exactly known results, or numerical results obtained by a finite-element method [68]. We also implemented the Skorokhod integral method (SIM) considered as a benchmark state-of-the-art Monte Carlo technique [55]. In Appendix D, we briefly describe the original method as being introduced by Schumm and Bressloff (and referred to as SIM1), and a minor variation of this method (referred to as SIM2). Throughout all cases, we fix $D = 1$. The discussion of cases shown in Fig. 3b, 3d, 3e is re-delegated to Appendix E.

## 4.1. Disk

For a disk of radius $R$ (Fig. 3a), the probability density $\Psi(\ell, p|\mathbf{x}_0)$ is known explicitly [23, 24]:

$$\Psi(\ell, p|\mathbf{x}_0) = \pi_0 \delta(\ell) + (1 - \pi_0) \mu_0^{(p)} e^{-\mu_0^{(p)} \ell}, \quad (47)$$
where $\alpha = \sqrt{p/D}$ and

$$ \mu_0^{(p)} = \alpha \frac{I_1(\alpha R)}{I_0(\alpha R)}, \qquad (48a) $$

$$ \pi_0 = 1 - \frac{I_0(\alpha r_0)}{I_0(\alpha R)}. \qquad (48b) $$

Once again, we have a mixture of the Dirac distribution and the exponential distribution. For this distribution, it is easy to compute all the moments of the $\ell_\delta$, in particular,

$$ \langle \ell_\delta \rangle = (1 - \pi_0)/\mu_0^{(p)}, \qquad \sigma^2 = \mathrm{Var}\{\ell_\delta\} = (1 - \pi_0^2)/(\mu_0^{(p)})^2. \qquad (49) $$

We fix $p = 1$ and the initial point $r_0 = 0.5$ inside the unit disk ($R = 1$). As the first test, we compute the empirical mean $\langle \ell_\delta \rangle_{\text{emp}}$ by using $N = 10^6$ particles and evaluate its relative error with respect to $\langle \ell_\delta \rangle$ by three methods (SIM, FLA, CLA) and three values of the layer width $\varepsilon$. Table 2 presents these results, as well as the CPU time for each simulation (all simulations were run on the same MacBook Pro with 1.4GHz quad-core Intel Core i5 processor and 8GB of 2133MHz LPDDR3 onboard memory). First of all, the three methods have similar CPU times for a given $\varepsilon$, the SIM being slightly faster. Importantly, the CPU time grows linearly with $1/\varepsilon$ for the three methods. This is in sharp contrast to the $\ln(1/\varepsilon)$ scaling of the CPU time for absorbing boundaries. Indeed, according to Eq. (3), the number of encounters with the boundary layer of width $\varepsilon$ behaves as $\ell_\delta/\varepsilon$, requiring on average $\langle \ell_\delta \rangle / \varepsilon$ escape events for each simulated trajectory. The expected linear scaling of the CPU time with $1/\varepsilon$ was the main motivation of our work, which aims at computing the boundary local time faster and more accurately by allowing larger $\varepsilon$.

Let us now compare the accuracy of the empirical mean computed by three methods. The relative error of the benchmark SIM1 changes from -7.43% to -0.91% and then to 1.80% as $\varepsilon$ goes down from $10^{-1}$ to $10^{-2}$ and then to $10^{-3}$. It seems that there is an “optimal” width layer $\varepsilon$ that minimizes the error, whereas further decreases of $\varepsilon$ is counter-productive. Note that a similar trend can be noticed in Fig. 3b of [55], when the relative error of the mean first-passage time is minimal at $\varepsilon \approx 0.01$. Even though such accuracy can be sufficient for some applications, it significantly exceeds the relative statistical error, estimated as $\frac{\sigma}{\langle \ell_\delta \rangle \sqrt{N}} \approx 0.1\%$. In other words, the relative error cannot be explained by statistical fluctuations (which are present in any finite sample) and thus originates from the employed approximations (e.g., finite $\varepsilon$). The accuracy of our modification (SIM2) is worse for this example.

The situation is very different for the FLA, for which the relative error drops from 4.8% at $\varepsilon = 10^{-1}$ to 0.1% at $\varepsilon = 10^{-3}$. In other words, we managed to reduce the relative error by factor $\sim 20$ and to achieve the level of statistical errors by using almost the same CPU time as for the SIM1. The accuracy is further improved by CLA. Even for $\varepsilon = 10^{-1}$, the relative error of the empirical mean is 0.5%, i.e., one achieves a two-fold improvement in accuracy by a 7 times faster algorithm, as compared to the SIM1 with $\varepsilon = 10^{-2}$. We conclude that the CLA dramatically outperforms two other methods.

The comparison between three methods carries on for the probability density $\Psi(\ell, p|x_0)$ of the $\ell_\delta$. Figure 4 shows the ratio between the empirical probability density obtained from the simulated sample of size $N = 10^6$, and the exact result in Eq. (47); in both cases, we exclude the trajectories that never touched the boundary, for which $\ell_\delta = 0$; in other words, we compare only regular parts of the probability density. The SIM1 yields a systematic bias in the probability density at large $\ell$, which is attenuated by decreasing $\varepsilon$ but still present even for $\varepsilon = 10^{-3}$ (panel (a)). This bias is also present for the FLA at $\varepsilon = 10^{-1}$ but it is already removed for $\varepsilon = 10^{-2}$ and $\varepsilon = 10^{-3}$. Expectedly, there is no such a bias for the CLA, even at $\varepsilon = 10^{-1}$. A larger scatter of points at large $\ell$ is a consequence of an insufficient statistics for rare trajectories with large $\ell_t$. Similar results were obtained for other simple domains such as a circular layer, a sphere, and a spherical shell, see Appendix E. We only mention that our modification SIM2 yields more accurate results in three dimensions than the original method SIM1 (see Appendix D for further discussion).

## 4.2. Random pack of disks in 2D

We proceed the validation by computing the statistics of the boundary local time $\ell_\delta$ of reflected Brownian motion outside a random pack of circular holed (disks) of the same radius $R_0$, confined by a large circle of radius $R$. The positions of these holes were chosen randomly. In the considered example shown in Fig. 3c, we set $R = 1$, and used 10 small holes of radius $R_0 = 0.1$. In this setting, the probability density of $\ell_\delta$ is still given by the spectral
<table><caption>Table 2: Comparison of SIM, FLA and CLA with $N = 10^6$ and different $\epsilon$ for the unit disk (Fig. 3a) with $p = 1$ and $r_0 = 0.5$. The mean and standard deviation of $\ell_\delta$ are $\langle\ell_\delta\rangle = 1.8817$, and $\sigma = 2.2113$, so that the relative statistical error is $\frac{\sigma}{\langle\ell_\delta\rangle\sqrt{N}} \approx 0.12\%$.</caption><thead><tr><th>Method</th><th>ε</th><th>CPU time (s)</th><th>⟨ℓ<sub>δ</sub>⟩<sub>emp</sub></th><th>Rel. Error</th><th>Method</th><th>ε</th><th>CPU time (s)</th><th>⟨ℓ<sub>δ</sub>⟩<sub>emp</sub></th><th>Rel. Error</th></tr></thead><tbody><tr><td rowspan="3">SIM1</td><td>10<sup>-1</sup></td><td>7.88 × 10<sup>0</sup></td><td>1.7420</td><td>-7.43%</td><td rowspan="3">SIM2</td><td>10<sup>-1</sup></td><td>7.37 × 10<sup>0</sup></td><td>1.7305</td><td>-8.04%</td></tr><tr><td>10<sup>-2</sup></td><td>6.97 × 10<sup>1</sup></td><td>1.8989</td><td>-0.91%</td><td>10<sup>-2</sup></td><td>5.62 × 10<sup>1</sup></td><td>1.7943</td><td>-4.65%</td></tr><tr><td>10<sup>-3</sup></td><td>6.86 × 10<sup>2</sup></td><td>1.9156</td><td>1.80%</td><td>10<sup>-3</sup></td><td>5.39 × 10<sup>2</sup></td><td>1.8006</td><td>-4.31%</td></tr><tr><td rowspan="3">FLA</td><td>10<sup>-1</sup></td><td>9.09 × 10<sup>0</sup></td><td>1.7906</td><td>-4.84%</td><td rowspan="3">CLA</td><td>10<sup>-1</sup></td><td>1.09 × 10<sup>1</sup></td><td>1.8912</td><td>0.50%</td></tr><tr><td>10<sup>-2</sup></td><td>7.87 × 10<sup>1</sup></td><td>1.8707</td><td>-0.59%</td><td>10<sup>-2</sup></td><td>9.56 × 10<sup>1</sup></td><td>1.8860</td><td>0.22%</td></tr><tr><td>10<sup>-3</sup></td><td>7.28 × 10<sup>2</sup></td><td>1.8798</td><td>-0.10%</td><td>10<sup>-3</sup></td><td>9.18 × 10<sup>2</sup></td><td>1.8808</td><td>-0.05%</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 4: Ratio $\Psi_{\text{emp}}(\ell, p|x_0)$ over $\Psi(\ell, p|x_0)$ as a function of $\ell/\langle\ell_\delta\rangle$ for the unit disk (Fig. 3a), with $p = 1$ and $r_0 = 0.5$. The empirical PDF $\Psi_{\text{emp}}$ was estimated from $N = 10^6$ simulated values of $\ell_\delta$ by (a) SIM1, (b) FLA, (c) CLA, with different choices of the boundary layer width: $\epsilon = 10^{-1}$ (triangles), $\epsilon = 10^{-2}$ (circles), and $\epsilon = 10^{-3}$ (squares). The trajectories that never touched the boundary and thus yielded $\ell_\delta = 0$, were excluded.</figcaption></figure>

<table><caption>Table 3: Comparison of SIM, FLA, CLA with $N = 10^6$ and $\epsilon = 10^{-3}$ for the unit disk perforated by 10 circular holes (Fig. 3c) with $p = 1$ and $R_0 = 0.1$. The “reference” mean value $\langle\ell_\sigma\rangle$ was computed from the spectral expansion (44) whose elements being obtained by a finite-element method [68]: $\langle\ell_\delta\rangle = 2.4783$ for $x_0 = (1, 0)$; $\langle\ell_\delta\rangle = 1.9400$ for $x_0 = (0, 0)$.</caption><thead><tr><th>x<sub>0</sub></th><th>Method</th><th>CPU time (s)</th><th>⟨ℓ<sub>δ</sub>⟩<sub>emp</sub></th><th>Rel. Error</th><th>x<sub>0</sub></th><th>Method</th><th>CPU time (s)</th><th>⟨ℓ<sub>δ</sub>⟩<sub>emp</sub></th><th>Rel. Error</th></tr></thead><tbody><tr><td rowspan="4">(+1, 0)</td><td>SIM1</td><td>4.47 × 10<sup>3</sup></td><td>2.5204</td><td>1.70%</td><td rowspan="4">(0, 0)</td><td>SIM1</td><td>4.24 × 10<sup>3</sup></td><td>1.9771</td><td>1.91%</td></tr><tr><td>SIM2</td><td>3.74 × 10<sup>3</sup></td><td>2.3648</td><td>-4.58%</td><td>SIM2</td><td>3.61 × 10<sup>3</sup></td><td>1.8522</td><td>-4.53%</td></tr><tr><td>FLA</td><td>4.18 × 10<sup>3</sup></td><td>2.4702</td><td>-0.33%</td><td>FLA</td><td>3.98 × 10<sup>3</sup></td><td>1.9319</td><td>-0.42%</td></tr><tr><td>CLA</td><td>4.51 × 10<sup>3</sup></td><td>2.4746</td><td>-0.15%</td><td>CLA</td><td>4.33 × 10<sup>3</sup></td><td>1.9402</td><td>0.01%</td></tr></tbody></table>

expansion (44); however, there is no explicit formula for the eigenvalues and eigenfunctions of the Steklov problem. We therefore compute them by a finite-element method [68]. The spectral expansion was truncated to 21 terms, and the maximal mesh size was 0.005 (244614 triangles). The accuracy of this computation was checked by changing the truncation order and the mesh size. In the following, we consider the probability density $\Psi(\ell, p|x_0)$ from Eq. (44) as a benchmark solution, to which Monte Carlo simulations of the boundary local time $\ell_\delta$ spent on the outer large circle of radius $R$ will be compared with. In turn, the encounters with the circular obstacles were ignored. This choice allows us to illustrate the flexibility of the proposed method, which can access the statistics of encounters with either the whole boundary, or its subsets.

Table 3 summarizes the CPU time, the empirical mean $\langle\ell_\delta\rangle_{\text{emp}}$, and its relative error for three methods. As the boundary layer width $\epsilon$ should necessarily be (much) smaller than the radius $R_0 = 0.1$ of circular holes, we provide results only for $\epsilon = 10^{-3}$. The conclusions on the comparison between three methods are quite similar to those for the unit disk: while the CPUs of three methods are comparable, the FLA and especially CLA provide much more accurate results. Figure 5 shows the comparison of empirical densities with those derived by the finite-element method, for $\epsilon = 10^{-3}$. As previously, the SIM1 yields a bias (which is consistent with the 2% relative error in the mean), whereas both FLA and CLA are not biased.
<figure><img src="image_6.png" /><figcaption>Figure 5: Ratio $\Psi_{\text{emp}}(\ell, p|x_0)$ over $\Psi(\ell, p|x_0)$ (obtained from the finite-element method) as a function of $\ell/\langle\ell_\delta\rangle$ for the unit disk perforated by 10 circular holes (Fig. 3c), with $p = 1$, $R_0 = 0.1$, and the boundary layer width $\varepsilon = 10^{-3}$. The empirical PDF $\Psi_{\text{emp}}$ was estimated from $N = 10^6$ simulated values of $\ell_\delta$ for different starting points: (a) $x_0 = (1, 0)$, (b) $x_0 = (0, 0)$, with different Monte Carlo methods: SIM1 (triangles), FLA (circles), and CLA (squares). The trajectories that never touched the boundary and thus yielded $\ell_\delta = 0$, were excluded.</figcaption></figure>

<table><caption>Table 4: Comparison between SIM and FLA with $N = 10^6$ and different $\varepsilon$ for the torus (Fig. 3f) with $p = 1$, $R_0 = 0.5$, $R = 1.5$, and the initial position (1.5, 0, 0). The “reference” mean value $\langle\ell_\delta\rangle = 3.8748$ was computed from the spectral expansion (44) whose elements being obtained by a finite-element method [68].</caption><thead><tr><th>Method</th><th>ε</th><th>CPU time (s)</th><th>⟨ℓ<sub>δ</sub>⟩<sub>emp</sub></th><th>Rel. Error</th><th>Method</th><th>ε</th><th>CPU time (s)</th><th>⟨ℓ<sub>δ</sub>⟩<sub>emp</sub></th><th>Rel. Error</th></tr></thead><tbody><tr><td rowspan="3">SIM1</td><td>10<sup>-1</sup></td><td>1.49 × 10<sup>1</sup></td><td>3.3024</td><td>-14.77%</td><td rowspan="3">SIM2</td><td>10<sup>-1</sup></td><td>1.55 × 10<sup>1</sup></td><td>3.5039</td><td>-9.57%</td></tr><tr><td>10<sup>-2</sup></td><td>2.38 × 10<sup>2</sup></td><td>4.0124</td><td>3.55%</td><td>10<sup>-2</sup></td><td>2.08 × 10<sup>2</sup></td><td>3.7825</td><td>-2.38%</td></tr><tr><td>10<sup>-3</sup></td><td>2.53 × 10<sup>3</sup></td><td>4.1076</td><td>6.01%</td><td>10<sup>-3</sup></td><td>2.06 × 10<sup>3</sup></td><td>3.8179</td><td>-1.47%</td></tr><tr><td rowspan="3">FLA</td><td>10<sup>-1</sup></td><td>2.86 × 10<sup>1</sup></td><td>3.2777</td><td>-15.41%</td><td rowspan="3"></td><td></td><td></td><td></td><td></td></tr><tr><td>10<sup>-2</sup></td><td>3.74 × 10<sup>2</sup></td><td>3.8034</td><td>-1.84%</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>10<sup>-3</sup></td><td>3.81 × 10<sup>3</sup></td><td>3.8629</td><td>-0.31%</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

### 4.3. Torus

As the last example, we consider a torus in three dimensions, which can be parameterized as:

$$ \begin{cases} x(\theta, \varphi) = (R + R_0 \cos \theta) \cos \varphi, \\ y(\theta, \varphi) = (R + R_0 \cos \theta) \sin \varphi, \\ z(\theta, \varphi) = R_0 \sin \theta, \end{cases} \qquad (50) $$

using angular coordinates $\theta, \varphi \in [0, 2\pi)$, where the minor radius $R_0$ is the radius of the tube, and the major radius $R$ is the distance between the center of the tube and the center of the torus. We select the torus with $R = 1.5$ and $R_0 = 0.5$ (Fig. 3f), and set $p = 1$. For this domain, we compute the eigenvalues and eigenfunctions of the Steklov problem by a finite-element method (adapted for axisymmetric three-dimensional problems) [68]. The spectral expansion (44) is truncated to 4 terms, and the maximal mesh size was 0.005 (75314 triangles). This benchmark solution is compared to the empirical probability density obtained by Monte Carlo simulations.

As the boundary points on the torus have two different radii of curvature, the CLA is not applicable. We test therefore only the FLA and the SIM. Table 4 compares the empirical means with the “reference” one. In this three-dimensional setting, our modification SIM2 yields more accurate results than the original method SIM1. Curiously, the FLA at $\varepsilon = 10^{-1}$ yields a bigger relative error than the SIM2. However, as $\varepsilon$ decreases, the relative error of the FLA drops much faster than that of both SIM1 and SIM2. Figure 6 presents the comparison of the PDFs of the boundary local time $\ell_\delta$, with similar conclusions.
<figure><img src="image_7.png" /><figcaption>Figure 6: Ratio $\Psi_{\text{emp}}(\ell, p|x_0)$ over $\Psi(\ell, p|x_0)$ (obtained by the finite-element method) as a function of $\ell/\langle\ell_\delta\rangle$ for the torus (Fig. 3f), with $p=1$, $R_0 = 0.5$, $R=1.5$, and the initial point located at (1.5, 0, 0). The empirical PDF $\Psi_{\text{emp}}$ was estimated from $N = 10^6$ simulated values of $\ell_\delta$ by (a) SIM1, (b) FLA, with different choices of the boundary layer width: $\varepsilon = 10^{-1}$ (triangles), $\varepsilon = 10^{-2}$ (circles), and $\varepsilon = 10^{-3}$ (squares). The trajectories that never touched the boundary and thus yielded $\ell_\delta = 0$, were excluded.</figcaption></figure>

## 5. Discussion and conclusion

In this paper, we proposed an efficient numerical tool to access the statistics of the boundary local time $\ell_t$ of reflected Brownian motion in Euclidean domains with smooth boundaries. Using the approximate solution of the escape problem for a thin boundary layer, we managed to replace a refined simulation of the random trajectory near the boundary by a single escape from that layer. This step eliminated the most time-consuming part of former Monte Carlo techniques and allowed us to relax the requirement of a too small layer width $\varepsilon$ that was needed to control the accuracy of the estimated boundary local time. Our implementation of this step employed several approximations: (i) the boundary layer was locally approximated by either a flat layer (FLA), or a curved layer with a known curvature (CLA); (ii) the conditional escape times $\tau_0$ and $\tau_\ell$ were replaced by their mean values $\langle\tau_0\rangle$ and $\langle\tau_\ell\rangle$; (iii) the escape position $X'_\tau$ was replaced by its mean value $\langle X'_\tau\rangle$; (iv) the exact stopping condition (e.g., $t_{k+1} = T$) was replaced by a loose condition (e.g., $t_{k+1} > T$). Despite these approximations, both FLA and CLA provided more accurate results than the state-of-the-art SIM under a comparable computational time. Moreover, the CLA was remarkably accurate even at relatively large boundary layer width (e.g., $\varepsilon/R = 0.1$). Even though the approximations (ii)-(iv) can be further relaxed (see discussions in Sec. 2.5 and Appendix B), the accuracy of the CLA was high enough to ignore these improvements.

At the same time, the approximation (i) was central to our approach. While the FLA is valid for any Euclidean domain in $\mathbb{R}^d$ ($d \ge 2$) with a smooth boundary, the applicability of the CLA is more restricted. In fact, it is valid for any planar domain with a smooth boundary, as well as for three-dimensional domains whose smooth boundary can be locally approximated by a sphere (i.e., with two equal radii of curvature). For instance, the CLA was not applicable for a torus that we discussed in Sec. 4.3. Further investigations of the escape problem from ellipsoidal or even more sophisticated layers, as well as to higher dimensions, present an interesting perspective of this work. Another important research direction is an extension of our approach to domains with angles (like polygons or prefractal domains) or more sophisticated singularities (e.g., cusps). In this case, one can consider escape events from triangular or rectangular regions, in analogy to efficient extensions of the WOS algorithm [39, 42]. For this purpose, one would need to find an approximate solution of the underlying escape problem [58].

As explained in Sec. 1, the boundary local time plays the central role in the theory of stochastic processes and has numerous applications. First of all, our numerical approach allows one to simulate reflected Brownian motion $X_t$ and the associated boundary local time $\ell_t$ in various confinements. For instance, random packs of monodisperse or polydisperse disks or spheres (as an example considered in Sec. 4.2) are typical models of porous media in material sciences, chemistry, and biology. In this way, one can compute the statistics of encounters between the diffusing particle and the boundary or a specific region on that boundary (e.g., we considered separately the boundary local times on the inner and outer circles of a circular annulus in Appendix E.1). Moreover, one can easily treat
more general problems with multiple diffusing particles or/and multiple targets, and access the joint distribution of various boundary local times and their mutual correlations [14, 69]. These quantities help to characterize the competition between particles for a target or the competition between targets for a particle. An extension of this numerical approach for simulating different functionals of the boundary local time presents an interesting perspective [7, 70]. From a broader perspective, one can study the boundary local time of a general diffusive process with given drift and covariance matrix [2, 3]. There are numerous Monte Carlo methods for simulating such stochastic processes as Markov chains at small timesteps (see [71–77] and references therein). While these methods can deal with much more general diffusive processes, their implementation in multi-scale media may be computationally prohibited due to small timesteps. However, if an applied potential and the induced drift are localized near the boundary, a combination of the above methods with the WOS algorithm can be beneficial. In this way, one could investigate the effect of the boundary potential (e.g., short-ranged electric attraction or repulsion) onto the boundary local time and the consequent surface reactions.

While the statistics of the boundary local time can be accessed by different numerical methods (e.g., by a finite-element method used in Sec. 4.2), Monte Carlo simulations have numerous advantages:

(i) One can directly generate the boundary local time $\ell_t$ at a fixed time $t$, whereas the methods relying on the spectral expansion (44) are less suitable in time domain, as they would require the inverse Laplace transform. For this reason, we did not look at the probability density function of $\ell_t$ that will be reported elsewhere.

(ii) Monte Carlo techniques are very flexible and can be easily implemented in any dimension; in contrast, an implementation of a finite-element method in dimensions $d \ge 3$ can be rather sophisticated and very time-consuming.

(iii) One can easily treat multiple targets (as the example shown in Fig. 3c) and compute distributions of the boundary local times $\ell_t^i$ on these targets. Such joint distributions are not easily accessible via other methods (see further discussions on multiple targets in [69]).

Potential applications of the proposed method go far beyond the statistics of encounters. For instance, the position $X_T$ of the process at the stopping time $\mathcal{T} = \inf\{t > 0 : \ell_t > \hat{\ell}\}$ when the boundary local time $\ell_t$ exceeds a random threshold $\hat{\ell}$, determines the position of the reaction event. The exponential distribution of the threshold, $\mathbb{P}\{\hat{\ell} > \ell\} = e^{-q\ell}$, corresponds to the most common case of a constant partial reactivity $\kappa = qD$ in diffusion-controlled reactions [15]. In this case, the distribution of the position $X_T$ is known as the spread harmonic measure [37, 78], whereas the stopping time $\mathcal{T}$ is called the first-reaction time or the first-passage time to the reaction event [79–85]. Both quantities can be accessed by the proposed method. However, other choices for the threshold distribution allow one to represent more sophisticated surface reactions such as progressive passivation or activation of the reactive boundary by the diffusing particles [4]. Our numerical method allows one to investigate the distribution of the reaction positions and thus to identify the most prolific subsets of the boundary with the highest reaction probability. This information can be further used to address various optimization problems in chemical engineering. From an even more general perspective, as Monte Carlo techniques are often employed to solve partial differential equations, the escape-from-a-layer approach can thus be helpful in treating problems with Neumann and Robin boundary conditions [25–27, 54, 66, 86–92].

## Acknowledgements

D.S.G. acknowledges a partial financial support from the Alexander von Humboldt Foundation through a Bessel Research Award.

## Appendix A. Escape time distribution

The WOS algorithm required generation of the random escape time $\tau_c$ from a disk or a ball. As this is a fairly standard procedure (see [93] and references therein), we only sketch the main steps. The survival probability of a particle located in the center of the unit disk, $S_c$, reads

$$S_c(t) = \mathbb{P}\{\tau_c > t\} = 2 \sum_{k=0}^{\infty} \frac{\exp(-\alpha_{0k}^2 t)}{\alpha_{0k} J_1(\alpha_{0k})}, \quad (A.1)$$
where $\alpha_{0k}$ refers to the k-th zero of the Bessel function $J_0(x)$ of the first kind, and we set $D = 1$. In turn, for the unit sphere with a particle located at center, the survival probability is

$$S_c(t) = 2 \sum_{k=1}^{\infty} (-1)^{k+1} e^{-\pi^2 k^2 t}. \qquad (A.2)$$

Explicit formulas are also known in higher dimensions.

In order to generate the random variable $\tau_c$, one needs to invert the survival probability, i.e., to find the function $T_c(x)$ such that $S_c[T_c(x)] = x$, for any $x \in (0, 1)$. If $\eta$ is uniformly distributed on $(0, 1)$, then $\tau_c = T_c(\eta)$ is the escape time. In practice, one can discretize the interval $(0, 1)$ into a sequence of points $x_k = k/K$, such as $k = 1, 2, \dots, K-1$ with large enough $K$, solve the equation $S_c(t_k) = x_k$ for each $k$ (e.g. by the bisection method), and set $T_c(x_k) = t_k$. Since the survival probability is a smooth, monotonously decreasing function, this computation can be done rapidly even for large $K$. Most importantly, this computation has to be done only once and then the stored values of $T_c(x_k)$ can be preloaded before running Monte Carlo simulations. Once a uniform random variable $\eta \in (0, 1)$ is generated, one rapidly gets $\tau_c = T_c(\eta)$ by interpolation.

In order to be able to accurately generate large escape times, one can rely on the long-time behavior:

$$S_c(t) \approx \begin{cases} \frac{2}{\alpha_{00} J_1(\alpha_{00})} e^{-\alpha_{00}^2 t} & (d=2), \\ 2e^{-\pi^2 t} & (d=3). \end{cases} \qquad (A.3)$$

If the random variable $\eta$ is below $x_1 = 1/K$ (i.e., $\tau$ is large), one can use Eq. (A.3) to invert $S_c(t)$ at large $t$, yielding an extrapolation relation

$$\tau_c = \begin{cases} -\frac{1}{\alpha_{00}^2} \ln \left[ \frac{\eta \alpha_{00} J_1(\alpha_{00})}{2} \right] & (d=2), \\ -\frac{1}{\pi^2} \ln \left( \frac{\eta}{2} \right) & (d=3). \end{cases} \qquad (A.4)$$

Replacing $\tau_c$ by $\rho^2/\tau_c$ allows one to get the escape time from a disk or a ball of radius $\rho$ and a given diffusion coefficient $D$.

To skip the above inversion of the survival probability, it is quite common to replace the random escape time $\tau_c$ by its mean value, which is particularly simple for the escape from the center of a ball of radius $\rho$ in $\mathbb{R}^d$: $\langle \tau_c \rangle = \frac{\rho^2}{2dD}$. However, this simplification can lead to significant errors when the particle starts far from the boundary so that the first jump to a large distance $\rho$ gives a large (fixed) contribution $\rho^2/(2dD)$ to the time counter (see further discussion in Appendix C). For this reason, we generated random escape times $\tau_c$ in all simulations of this paper.

## Appendix B. Conditional probability density

In this Appendix, we briefly discuss the difficulties in generating conditional escape times and their possible solutions. We focus here on the escape from a flat layer and consider separately two conditional escape times: without hitting the bottom surface and with reflections on it.

The inverse Laplace transform of Eq. (13) yields

$$H_0(t|y_0) = \frac{2\pi D}{\varepsilon y_0} \sum_{k=1}^{\infty} (-1)^{k-1} k \sin(\pi k y_0 / \varepsilon) e^{-D \pi^2 k^2 / \varepsilon^2}. \qquad (B.1)$$

The survival probability of this random variable,

$$S_0(t|y_0) = \mathbb{P}_{y_0} \{\tau_0 > t\} = \int_t^\infty dt' H_0(t'|y_0) = \frac{2\varepsilon}{\pi y_0} \sum_{k=1}^\infty \frac{(-1)^{k-1}}{k} \sin(\pi k y_0 / \varepsilon) e^{-D \pi^2 k^2 / \varepsilon^2}, \qquad (B.2)$$
can be numerically inverted to generate the conditional escape time $\tau_0$ (see Appendix A for a similar procedure for a disk). However, the main practical difficulty is that this survival probability depends on $y_0/\epsilon$ so that its inversion may be quite costly.

The situation is much more difficult for the second conditional escape time $\tau_\ell$, whose moment-generating function is given by Eq. (15). There are two technical difficulties to perform its inverse Laplace transform: first, there is no simple explicit representation for this inversion (see, however, some related representations derived in [69]); second, there are three parameters, $t$, $y_0$ and $\ell$, that determine this probability density. As a consequence, a direct computation of this quantity for the whole set of parameters could be problematic.

To partly overcome the second difficulty, one can employ the following simplification. The expression in Eq. (15) can be written as

$$ \tilde{H}(p|y_0, \ell) = \frac{\epsilon \sinh(\sqrt{p/D}(\epsilon - y_0))}{(\epsilon - y_0) \sinh(\sqrt{p/D}\epsilon)} \tilde{H}(p|0, \ell). \quad (B.3) $$

The first factor is the Laplace transform of the probability density of the escape time, conditioned to the escape through the endpoint 0, which is actually equal to $H_0(t|\epsilon - y_0)$ due to the symmetry. In turn, the second factor is the Laplace transform of the conditional probability density $H(t|0, \ell)$ when the starting point is at 0. The product of these factors indicates that the escape time $\tau_\ell$ can be written as the sum of two independent random times: the time $\tau_1$ to reach the endpoint 0 and the time $\tau_2$ to cross the interval $(0, \epsilon)$ and to escape through the endpoint $\epsilon$. The probability density of the latter time is then determined by the inverse Laplace transform:

$$ H(t|0, \ell) = \mathcal{L}^{-1}\{\tilde{H}(p|0, \ell)\} = \mathcal{L}^{-1}\left\{\frac{\epsilon \sqrt{p/D}}{\sinh(\sqrt{p/D}\epsilon)} e^{-\ell(\sqrt{p/D}\tanh(\sqrt{p/D}\epsilon) - 1/\epsilon)}\right\}. \quad (B.4) $$

Even though this inverse Laplace transform needs to be computed numerically, it does not depend on $y_0$ that simplifies its practical implementation. In the same way, one can compute the related survival probability:

$$ S(t|0, \ell) = 1 - \mathcal{L}^{-1}\left\{\frac{\tilde{H}_0(p|0, \ell)}{p}\right\} = 1 - \mathcal{L}^{-1}\left\{\frac{\epsilon \sqrt{p/D}}{p \sinh(\sqrt{p/D}\epsilon)} e^{-\ell(\sqrt{p/D}\tanh(\sqrt{p/D}\epsilon) - 1/\epsilon)}\right\}. \quad (B.5) $$

To avoid technical issues with Laplace transform inversion, we did not use these exact representations and replaced the random variables $\tau_0$ and $\tau_\ell$ by their mean values, which were computed explicitly.

## Appendix C. Justification of the approximation by means

One practically relevant approximation consists in replacing random escape times by their mean values. At a first thought, this approximation may sound rather crude. While a rigorous analysis of the accuracy of this approximation is beyond the scope of the manuscript, we provide some qualitative arguments and rationales.

To convey the idea, we consider a simpler, slightly different problem of estimating the first-passage time $T$ to an absorbing boundary and assume that a random trajectory is simulated by generating a sequence of spheres of the same radius $\rho$. In other words, one draws a sphere of radius $\rho$ centered at the starting position and moves the particle to a random point on that sphere. The escape time from the sphere, $\tau$, is a random variable, with explicitly known distribution (see Appendix A). The process is repeated again and again, until the newly drawn $(n+1)\text{th}$ sphere crosses the boundary; the number $n$ of generated spheres is thus also a random variable. The first-passage time $T$ to the boundary is then estimated by summing independent random escape times generated until hitting the boundary: $T = \tau_1 + \tau_2 + \dots + \tau_n$. In turn, the approximation consists in replacing random variables $\tau_1, \tau_2, \dots, \tau_n$ by their means, $\langle \tau \rangle = \rho^2/(6D)$, i.e., $T_{\text{app}} = n\langle \tau \rangle = n\rho^2/(6D)$. While this approximation yields the correct mean FPT, $\langle T \rangle = \langle n \rangle \langle \tau \rangle = \langle T_{\text{app}} \rangle$, the distributions of $T$ and $T_{\text{app}}$ are different. How accurately can the distribution of $T_{\text{app}}$ approximate that of $T$?

To estimate the error of this approximation, let us compute the variance of $T$. A straightforward computation gives:

$$ \mathrm{Var}\{T\} = \langle n \rangle \mathrm{Var}\{\tau\} + \mathrm{Var}\{n\}\langle \tau \rangle^2. \quad (C.1) $$
<figure><img src="image_8.png" /><figcaption>Figure D.7: Illustration of the reflection on a circular boundary. When a particle (red dot) is inside the boundary layer (light orange) of width $\epsilon$, one draws a circle of radius $\rho = 2\epsilon$ and generates a random position on it. If this position (green dot) turns out to be outside the domain, the final position (blue dot) is either attached on the boundary along the normal direction (violet dashed line) in SIM1 (a), or obtained by a mirror reflection along the normal direction (violet dashed line) in SIM2 (b).</figcaption></figure>

As the approximation replaces random variables $\tau_1, \tau_2, \dots, \tau_n$ by their means, the variance of the approximated first-passage time is $\text{Var}\{T_{\text{app}}\} = \text{Var}\{n\}\langle\tau\rangle^2$, i.e., it ignores the variance of $\tau$ (the first term in Eq. (C.1)). We conclude that the proposed approximation always underestimates the variance of the first-passage time.

However, in many practical cases, the second term in Eq. (C.1) is dominant so that $\text{Var}\{T\} \approx \text{Var}\{T_{\text{app}}\}$, suggesting that the approximation can be accurate. Indeed, the variance of the escape time from a sphere, $\text{Var}\{\tau\}$, can be found explicitly; but even without an exact computation, it must be proportional to $(\rho^2/D)^2$ from dimensional analysis, and thus to $\langle\tau\rangle^2$. Moreover, the variance of the FPT is also generally proportional to the square of its mean, $\text{Var}\{T\} \sim \langle T\rangle^2$, i.e., Eq. (C.1) reads as $\langle n\rangle^2 \sim \langle n\rangle + \text{Var}\{n\}$. If the mean number of generated spheres, $\langle n\rangle$, is large enough, the variance $\text{Var}\{n\}$ scales as $\langle n\rangle^2$, i.e., the relative contribution of the first term to Eq. (C.1) is of the order of $1/\langle n\rangle$ and thus can be neglected.

We stress that the above arguments are qualitative and may not be valid in general. For instance, the relation $\text{Var}\{T\} \sim \langle T\rangle^2$ may not hold. Moreover, when the radii of generated spheres are not fixed (as in the WOS algorithm), the contribution from few large spheres can be dominant. For this reason, we do not apply the approximation to simulate Brownian motion in the bulk, i.e., when using spheres of variable radii; in this part of the algorithm, we generate random escape times (see Appendix A). The approximation is only used to simulate the escape from a layer of a fixed small width $\epsilon$.

## Appendix D. Skorokhod integral method

For completeness, we briefly recall here the Skorokhod integral method developed in [55]. The WOS algorithm is employed until the particle enters the boundary layer of width $\epsilon$. Then, the radius $\rho$ is fixed to be proportional to $\epsilon$: $\rho = 2\epsilon$ in two dimensions, and $\rho = 3\epsilon$ in three dimensions. After the jump to the distance $\rho$, the particle can leave the boundary layer, stay inside it, or exit the domain. In the first case, the WOS algorithm with variable-distance jumps is resumed; in the second case, one makes another jumps with the radius $\rho = d\epsilon$. In turn, the latter option is treated as a reflection event. The current position generated outside the domain is reflected back with respect to the boundary (Fig. D.7a). The time counter is increased by $\delta = \frac{\ell_D^2}{D}\tau_c$, where $\tau_c$ is the escape time from the unit disk or the unit ball. In turn, the Skorokhod integral representation allows one to increment the boundary local time as [54, 55]:

$$ \ell_{\ell_{k+1}} = \ell_{\ell_k} + \sqrt{\pi/2} \sqrt{D\delta}. \qquad (D.1) $$

We refer to this method as SIM1.

Despite a number of tests realized in [54, 55], we provide an independent validation to check that our implementation of the SIM1 yields the correct results. In particular, we check and validate the choice $\rho = d\epsilon$ for the radius $\rho$ inside the boundary layer. For this purpose, we compute numerically the spread harmonic measure density $\omega_q(x|x_0)$ on the circle and on the sphere, for which exact solutions are available [78]. For a fixed starting point $x_0$, $\omega_q(x|x_0)d\mathbf{x}$ is the probability that $X_T \in (x, x + d\mathbf{x})$ where $T = \inf\{t > 0 : \ell_t > \hat{\ell}\}$ is the first crossing time of a random threshold $\hat{\ell}$ by the $\ell_t$. In other words, we employ here the third stopping condition (see Sec. 2.5) with the exponential threshold $\hat{\ell}$ such that $P\{\hat{\ell} > \ell\} = e^{-q\ell}$, for a given $q > 0$. From the applicative point of view, $X_T$ is the random position of the reaction event on the boundary with partial reactivity, $\kappa = qD$ [4, 24].
<figure><img src="image_9.png" /><figcaption>Figure D.8: (a) Spread harmonic measure density on the unit circle ($R = 1$): comparison between the analytical expression (D.3) with $K = 60$, truncated at $j_{max} = 50$, shown by the black line, and numerical estimates by SIM1. Parameters are: the boundary layer width $\epsilon = 10^{-3}$, reactivity $q = 10$, particle number $N = 10^6$, and the initial point is located at (0, 0.9), with different choices of the ratio $\rho/\epsilon$: $\rho = \epsilon$ (triangles), $\rho = 2\epsilon$ (squares), and $\rho = 3\epsilon$ (circles). (b) Spread harmonic measure density on the unit sphere ($R = 1$): comparison between the analytical expression (D.6) with $K = 60$, truncated at $n_{max} = 60$, shown by black line, and numerical ones by SIM1. Parameters are: the boundary layer width $\epsilon = 10^{-3}$, reactivity $q = 10$, particle number $N = 10^7$, and the initial point is located at (0, 0, 0.9), with different choices of the ratio $\rho/\epsilon$: $\rho = 2\epsilon$ (squares), $\rho = 3\epsilon$ (circles), and $\rho = 4\epsilon$ (triangles). (c) Ratio of numerical estimates $p_k$ shown in panel (a) to the analytical results (D.3), with different choices of the ratio $\rho/\epsilon$: $\rho = \epsilon$ (triangles), $\rho = 2\epsilon$ (squares), and $\rho = 3\epsilon$ (circles). (d) Ratio of numerical estimates $p_k$ shown in panel (b) to the analytical results (D.6), with different choices of the ratio $\rho/\epsilon$: $\rho = 2\epsilon$ (squares), $\rho = 3\epsilon$ (circles), and $\rho = 4\epsilon$ (triangles). (e) Ratio of numerical estimates $p_k$ obtained by SIM2 (not shown) to the analytical results (D.3), with different choices of the ratio $\rho/\epsilon$: $\rho = \epsilon$ (triangles), $\rho = 2\epsilon$ (squares), and $\rho = 3\epsilon$ (circles). (f) Ratio of numerical estimates $p_k$ obtained by SIM2 (not shown) to the analytical results (D.6), with different choices of the ratio $\rho/\epsilon$: $\rho = 2\epsilon$ (squares), $\rho = 3\epsilon$ (circles), and $\rho = 4\epsilon$ (triangles).</figcaption></figure>

## Appendix D.1. Test for a disk

The spread harmonic measure density on the circle of radius $R$ reads [78]

$$ \omega_q(\theta|r_0, \theta_0) = \frac{1}{2\pi R} \left( 1 + 2 \sum_{j=1}^{\infty} \left(\frac{r_0}{R}\right)^j \frac{\cos[j(\theta - \theta_0)]}{1 + \frac{j}{qR}} \right), \qquad (D.2) $$

where $(r_0, \theta_0)$ is the starting point in polar coordinates. In order to estimate this density from Monte Carlo simulations, we divide the circle into $K$ equal arcs such that the $k$-th arc contains $\theta \in [k-1, k) \times \frac{2\pi}{K}$, $k = 1, 2, 3, \dots, K$. The probability $p_k$ of the reaction event on this arc reads then

$$ p_k = R \int_{2\pi(k-1)/K}^{2\pi k/K} \omega_q(\theta|r_0, \theta_0) d\theta = \frac{1}{K} + \sum_{j=1}^{\infty} \frac{2qR}{\pi j(j+qR)} \left(\frac{r_0}{R}\right)^j \sin\left(\frac{\pi j}{K}\right) \cos\left(j\left[\frac{2\pi}{K}\left(k-\frac{1}{2}\right) - \theta_0\right]\right). \qquad (D.3) $$

In turn, the probabilities $p_k$ can be estimated by running Monte Carlo simulations with the third stopping condition and counting the number $N_k$ of trajectories that stopped on the $k$-th arc, so that $N_k/N$ approximates $p_k$.

Figure D.8a compares the exact probabilities $p_k$ from Eq. (D.3) to the empirical ones obtained by Monte Carlo simulations with SIM1. We compare three choices for the radius $\rho$ inside the boundary layer and illustrate that $\rho = 2\epsilon$ yields the most accurate results. Other starting points and values of $q$ were tested and gave similar results (not shown).
## Appendix D.2. Test for a sphere

For a sphere of radius $R$, the spread harmonic measure density reads [94]:

$$ \omega_q(\mathbf{x}|\mathbf{x}_0) = \frac{1}{4\pi R^2} \sum_{n=0}^{\infty} P_n(\cos \gamma) \left(\frac{|\mathbf{x}_0|}{R}\right)^n \frac{2n+1}{1+n/(qR)}, \qquad (D.4) $$

where $\mathbf{x}$ is the position on the sphere, $\cos \gamma = \frac{\langle \mathbf{x} \cdot \mathbf{x}_0 \rangle}{|\mathbf{x}| |\mathbf{x}_0|}$, $P_n(\mathbf{x})$ is the Legendre polynomial. Using the spherical coordinates, $\mathbf{x} = (R, \theta, \phi)$, $\mathbf{x}_0 = (|\mathbf{x}_0|, \theta_0, \phi_0)$, it is convenient to integrate over the angle $\phi$ to get

$$ \bar{\omega}(\theta|\theta_0) = \frac{1}{2} \sum_{n=0}^{\infty} P_n(\cos \theta) P_n(\cos \theta_0) \left(\frac{r_0}{R}\right)^n \frac{2n+1}{1+n/(qR)}. \qquad (D.5) $$

As previously, we divide the sphere into $K$ regions such that the $k$-th layer contains $\theta \in [k-1, k) \times \frac{\pi}{K}$, with $k = 1, 2, \dots, K$. The probability $p_k$ of the reaction event on $k$-th region is thus expressed by

$$ p_k = \int_{(k-1)\pi/K}^{k\pi/K} \bar{\omega}(\theta|\theta_0) \sin\theta d\theta = \frac{\xi_{k-1} - \xi_k}{2} - \frac{1}{2} \sum_{n=1}^{\infty} \frac{(r_0/R)^n}{1+n/(qR)} \{[P_{n+1}(\xi_k) - P_{n-1}(\xi_k)] - [P_{n+1}(\xi_{k-1}) - P_{n-1}(\xi_{k-1})]\}, \qquad (D.6) $$

with $\xi_k = \cos(k\pi/K)$.

We test different choices for the ratio $\rho/\epsilon$ by computing the spread harmonic measure on a unit sphere with $R = 1, \epsilon = 10^{-3}, q = 10$, and the initial point is located at $(0, 0, 0.9)$. Figure D.8b compares the exact probabilities $p_k$ from Eq. (D.6) to the empirical ones obtained by Monte Carlo simulations with SIM1. In particular, panel (d) indicates that the choice $\rho = 3\epsilon$ suggested in [43], does not yield the most accurate results. Moreover, the results remain biased even for a larger value $\rho = 4\epsilon$. This bias can potentially explain higher relative errors of the SIM1 in three dimensions.

## Appendix D.3. A modification of SIM

While Schumm and Bressloff designed and validated their SIM for planar domains [55], its extension to three-dimensional domains was straightforward. However, no validation of this extension was performed. Our test for a sphere in Appendix D.2 revealed the accuracy issues of this extension. Moreover, the relative error of the mean $\langle l_\delta \rangle$ was quite big for all the considered three-dimensional domains, even at $\epsilon = 10^{-3}$ (see Tables 4, E.6, E.7). In an attempt to understand this discrepancy, we tried different ways of modeling the reflection event. In fact, the original method suggests to relocate the particle, which left the domain, to the boundary point along the normal direction (Fig. D.7a). To some extent, this is an arbitrary choice. For instance, most algorithms simulating reflected diffusion employ a mirror reflection across the boundary (see Fig. D.7b), but other choices are as well possible. It is generally believed that the choice of reflection, occurring at a small scale $\epsilon$, provides only a minor effect onto simulation results. However, as the boundary local time is incremented only at such reflections, the minor errors potentially induced by the “wrong” choice, can accumulate and result in considerable bias.

In order to check this point, we propose a minor modification of the SIM, which consists in two steps: (i) the relocation to the boundary is replaced by a mirror reflection (compare Fig. D.7a and D.7b); (ii) after the reflection, we accumulated the durations of the consecutive jumps inside the boundary layer, until the particle leaves it. The accumulated duration $\delta$ is then used to increment the boundary local time via Eq. (D.1). We refer to this method as SIM2. Figure D.8f shows that the SIM2 with $\rho = d\epsilon$ yields unbiased results for the spread harmonic measure in both two- and three-dimensional cases. We used this choice throughout all the manuscript. Tables 4, E.6, E.7 indicate that SIM2 yields more accurate results than SIM1 in three-dimensional examples, but less accurate results in two-dimensional examples. A more systematic analysis of this observation presents an interesting problem, which is, however, beyond the scope of this paper.
<table><caption>Table E.5: Comparison of the SIM, FLA, CLA with $N = 10^6$ and different $\epsilon$ for the annulus (Fig. 3b) with $p = 1$, $R = 1$, $L = 2$, and $r_0 = 1.5$. The mean and standard deviation are: $\langle \ell_{\delta}^{i} \rangle = 0.6494$, $\sigma^i = 0.9112$, and $\langle \ell_{\delta}^{o} \rangle = 1.2660$, $\sigma^o = 1.5931$, so that the relative statistical errors are $\frac{\sigma^i}{\langle \ell_{\delta}^{i} \rangle \sqrt{N}} \approx 0.14\%$, and $\frac{\sigma^o}{\langle \ell_{\delta}^{o} \rangle \sqrt{N}} \approx 0.13\%$.</caption><thead><tr><th>Method</th><th>ε</th><th>CPU time (s)</th><th>⟨ℓ̂<sup>i</sup><sub>δ</sub>⟩<sub>emp</sub></th><th>Rel. Error</th><th>⟨ℓ̂<sup>o</sup><sub>δ</sub>⟩<sub>emp</sub></th><th>Rel. Error</th></tr></thead><tbody><tr><td rowspan="3">SIM1</td><td>10<sup>-1</sup></td><td>9.57 × 10<sup>0</sup></td><td>0.5516</td><td>-15.07%</td><td>1.1524</td><td>-8.98%</td></tr><tr><td>10<sup>-2</sup></td><td>7.63 × 10<sup>1</sup></td><td>0.6504</td><td>0.15%</td><td>1.2769</td><td>0.86%</td></tr><tr><td>10<sup>-3</sup></td><td>7.11 × 10<sup>2</sup></td><td>0.6621</td><td>1.95%</td><td>1.2914</td><td>2.02%</td></tr><tr><td rowspan="3">SIM2</td><td>10<sup>-1</sup></td><td>8.87 × 10<sup>0</sup></td><td>0.6386</td><td>-1.66%</td><td>1.1820</td><td>-6.63%</td></tr><tr><td>10<sup>-2</sup></td><td>6.32 × 10<sup>1</sup></td><td>0.6217</td><td>-4.27%</td><td>1.2077</td><td>-4.61%</td></tr><tr><td>10<sup>-3</sup></td><td>5.73 × 10<sup>2</sup></td><td>0.6185</td><td>-4.75%</td><td>1.2091</td><td>-4.49%</td></tr><tr><td rowspan="3">FLA</td><td>10<sup>-1</sup></td><td>1.31 × 10<sup>1</sup></td><td>0.6883</td><td>6.00%</td><td>1.2344</td><td>-2.50%</td></tr><tr><td>10<sup>-2</sup></td><td>9.21 × 10<sup>1</sup></td><td>0.6532</td><td>0.58%</td><td>1.2619</td><td>-0.33%</td></tr><tr><td>10<sup>-3</sup></td><td>8.02 × 10<sup>2</sup></td><td>0.6508</td><td>0.22%</td><td>1.2665</td><td>-0.04%</td></tr><tr><td rowspan="3">CLA</td><td>10<sup>-1</sup></td><td>1.31 × 10<sup>1</sup></td><td>0.6522</td><td>0.43%</td><td>1.2694</td><td>0.27%</td></tr><tr><td>10<sup>-2</sup></td><td>1.04 × 10<sup>2</sup></td><td>0.6483</td><td>-0.17%</td><td>1.2645</td><td>-0.12%</td></tr><tr><td>10<sup>-3</sup></td><td>9.74 × 10<sup>2</sup></td><td>0.6483</td><td>-0.02%</td><td>1.2668</td><td>-0.06%</td></tr></tbody></table>

## Appendix E. Validation for other domains

### Appendix E.1. Annulus in 2D

We validate Monte Carlo simulations for a circular annulus of radii $R$ and $L$ (Fig. 3b). In this case, we evaluate separately the statistics of the boundary local times $\ell_{\delta}^{i}$ and $\ell_{\delta}^{o}$ on the inner and outer circles. The rotational symmetry of this domain implies again the simple form (47), with the following parameters [24]:

(i) for the boundary local time $\ell_{\delta}^{i}$ on the inner circle:

$$ \pi_0 = 1 - \frac{I_1(\alpha L)K_0(\alpha r_0) + K_1(\alpha L)I_0(\alpha r_0)}{I_1(\alpha L)K_0(\alpha R) + K_1(\alpha L)I_0(\alpha R)}, \qquad (\text{E.1a}) $$

$$ \mu_0^{(p)} = \alpha \frac{I_1(\alpha L)K_1(\alpha R) - K_1(\alpha L)I_1(\alpha R)}{I_1(\alpha L)K_0(\alpha R) + K_1(\alpha L)I_0(\alpha R)}. \qquad (\text{E.1b}) $$

(ii) for the boundary local time $\ell_{\delta}^{o}$ on the outer circle:

$$ \pi_0 = 1 - \frac{K_1(\alpha R)I_0(\alpha r_0) + I_1(\alpha R)K_0(\alpha r_0)}{K_1(\alpha R)I_0(\alpha L) + I_1(\alpha R)K_0(\alpha L)}, \qquad (\text{E.2a}) $$

$$ \mu_0^{(p)} = \alpha \frac{K_1(\alpha R)I_1(\alpha L) - I_1(\alpha R)K_1(\alpha L)}{K_1(\alpha R)I_0(\alpha L) + I_1(\alpha R)K_0(\alpha L)}. \qquad (\text{E.2b}) $$

As previously for the case of a disk, Table E.5 compares the empirical mean values of $\ell_{\delta}^{i}$ and $\ell_{\delta}^{o}$ found by different methods and different widths $\epsilon$ with their analytical expectations, given by Eq. (49). As in Table 2, the most accurate results by the SIM1 are achieved at an intermediate width $\epsilon = 10^{-2}$, whereas its accuracy does not improve at smaller $\epsilon$. For other methods, smaller $\epsilon$ results in more accurate values. The CLA outperforms two other methods for both $\ell_{\delta}^{i}$ and $\ell_{\delta}^{o}$. Accidentally, the relative error of $\langle \ell_{\delta}^{i} \rangle_{\text{emp}}$ for the SIM1 at $\epsilon = 10^{-2}$ is comparable to that of the CLA. However, this is not true for $\langle \ell_{\delta}^{o} \rangle_{\text{emp}}$.

Figure E.9 shows the comparison of the statistics of two boundary local times, $\ell_{\delta}^{i}$ and $\ell_{\delta}^{o}$, obtained by different Monte Carlo methods, with the exact solution (47), for $p = 1$ and several values of $\epsilon = 10^{-1}$, $10^{-2}$, $10^{-3}$. Curiously, the FLA presents opposite trends for $\ell_{\delta}^{i}$ and $\ell_{\delta}^{o}$ as $\ell$ increases. The CLA yields accurate results even for $\epsilon = 10^{-1}$.

### Appendix E.2. Sphere and spherical shell in 3D

The expression (47) is also valid for the sphere of radius $R$, with

$$ \mu_0^{(p)} = \alpha \operatorname{ctanh}(\alpha R) - \frac{1}{R}, \qquad \pi_0 = 1 - \frac{R \sinh(\alpha r_0)}{r_0 \sinh(\alpha R)}, \qquad (\text{E.3}) $$
<figure><img src="image_10.png" /><figcaption>Figure E.9: Ratio $\Psi_{\text{emp}}(\ell, p|x_0)$ over $\Psi(\ell, p|x_0)$ as a function of $\ell/\langle\ell_\delta\rangle$ for the circular annulus (Fig. 3b), with $p = 1$, $R = 1$, $L = 2$, and $r_0 = 1.5$. The empirical PDF $\Psi_{\text{emp}}$ was estimated from $N = 10^6$ simulated values of $\ell_\delta$ by (a,b) SIM1, (c,d) FLA, (e,f) CLA, with different choices of the boundary layer width: $\varepsilon = 10^{-1}$ (triangles), $\varepsilon = 10^{-2}$ (circles), and $\varepsilon = 10^{-3}$ (squares). (a,c,e) refer to $\ell_\delta^0$ on the inner circle and (b,d,f) refer to $\ell_\delta^1$ on the outer circle. The trajectories that never touched the boundary and thus yielded $\ell_\delta = 0$, were excluded.</figcaption></figure>

<table><caption>Table E.6: Comparison of the SIM, FLA, CLA with $N = 10^6$ and different $\varepsilon$ for the unit sphere (Fig. 3d) with $p = 1$ and $r_0 = 0.5$. The mean and standard deviation of $\ell_\delta$ are $\langle\ell_\delta\rangle = 2.8330$ and $\sigma = 3.1740$, so that the relative statistical error is $\frac{\sigma}{\langle\ell_\delta\rangle\sqrt{N}} \approx 0.11\%$.</caption><thead><tr><th>Method</th><th>ε</th><th>CPU time (s)</th><th>⟨ℓ<sub>δ</sub>⟩<sub>emp</sub></th><th>Rel. Error</th><th>Method</th><th>ε</th><th>CPU time (s)</th><th>⟨ℓ<sub>δ</sub>⟩<sub>emp</sub></th><th>Rel. Error</th></tr></thead><tbody><tr><td rowspan="3">SIM1</td><td>10<sup>-1</sup></td><td>8.23 × 10<sup>0</sup></td><td>2.5995</td><td>-8.24%</td><td rowspan="3">SIM2</td><td>10<sup>-1</sup></td><td>8.24 × 10<sup>0</sup></td><td>2.5394</td><td>-10.36%</td></tr><tr><td>10<sup>-2</sup></td><td>1.22 × 10<sup>2</sup></td><td>2.9557</td><td>4.33%</td><td>10<sup>-2</sup></td><td>1.04 × 10<sup>2</sup></td><td>2.7688</td><td>-2.27%</td></tr><tr><td>10<sup>-3</sup></td><td>1.33 × 10<sup>3</sup></td><td>3.0006</td><td>5.91%</td><td>10<sup>-3</sup></td><td>1.08 × 10<sup>3</sup></td><td>2.7893</td><td>-1.54%</td></tr><tr><td rowspan="3">FLA</td><td>10<sup>-1</sup></td><td>1.54 × 10<sup>1</sup></td><td>2.5353</td><td>-10.51%</td><td rowspan="3">CLA</td><td>10<sup>-1</sup></td><td>1.53 × 10<sup>1</sup></td><td>2.8476</td><td>0.52%</td></tr><tr><td>10<sup>-2</sup></td><td>1.82 × 10<sup>2</sup></td><td>2.7988</td><td>-1.21%</td><td>10<sup>-2</sup></td><td>1.90 × 10<sup>2</sup></td><td>2.8353</td><td>0.08%</td></tr><tr><td>10<sup>-3</sup></td><td>1.91 × 10<sup>3</sup></td><td>2.8351</td><td>0.07%</td><td>10<sup>-3</sup></td><td>2.01 × 10<sup>3</sup></td><td>2.8362</td><td>0.11%</td></tr></tbody></table>

where $r_0$ refers to the distance between the initial point and the sphere center, and $\alpha = \sqrt{p/D}$.

Figure E.10 presents the comparison of the exact PDF in Eq. (47) with the empirical ones obtained by different Monte Carlo methods. For $\varepsilon = 10^{-1}$, both SIM and FLA present a bias, which is removed at smaller $\varepsilon$. In turn, the CLA shows again excellent agreement for all cases.

Table E.6 compares three methods for evaluating the mean value $\langle\ell_\delta\rangle$. As for the disk, the SIM is outperformed by the FLA (especially for small $\varepsilon$), whereas the CLA yields the most accurate results. Curiously, the accuracies are higher in three dimensions than in two dimensions. Note that the SIM2 yields more accurate results than SIM1.

Furthermore, we also validate Monte Carlo simulations for a spherical shell between concentric spheres of radii $R$ and $L$ by evaluating $\ell_\delta^1$ and $\ell_\delta^0$. Their PDFs are still given by Eq. (47), with

(i) for the boundary local time $\ell_\delta^1$ on the inner sphere:

$$ \pi_0 = 1 - \frac{R \alpha L \cosh(\alpha(L - r_0)) - \sinh(\alpha(L - r_0))}{r_0 \alpha L \cosh(\alpha(L - R)) - \sinh(\alpha(L - R))}, \qquad (\text{E.4a}) $$

$$ \mu_0^{(p)} = \frac{(\alpha^2 LR - 1) \sinh(\alpha(L - R)) + \alpha(L - R) \cosh(\alpha(L - R))}{R(\alpha L \cosh(\alpha(L - R)) - \sinh(\alpha(L - R)))}. \qquad (\text{E.4b}) $$
<figure><img src="image_11.png" /><figcaption>Figure E.10: Ratio $\Psi_{\text{emp}}(\ell, p|x_0)$ over $\Psi(\ell, p|x_0)$ as a function of $\ell/\langle\ell_\delta\rangle$ for the unit sphere (Fig. 3d), with $p = 1$ and $r_0 = 0.5$. The empirical PDF $\Psi_{\text{emp}}$ was estimated from $N = 10^6$ simulated values of $\ell_\delta$ by (a) SIM1, (b) FLA, (c) CLA, with different choices of the boundary layer width: $\varepsilon = 10^{-1}$ (triangles), $\varepsilon = 10^{-2}$ (circles), and $\varepsilon = 10^{-3}$ (squares). The trajectories that never touched the boundary and thus yielded $\ell_\delta = 0$, were excluded.</figcaption></figure>

<table><caption>Table E.7: Comparison of the SIM, FLA, CLA with $N = 10^6$ and different $\varepsilon$ for the spherical shell (Fig. 3e) with $p = 1$, $R = 1$, $L = 2$, and $r_0 = 1.5$. The mean and standard deviation are: $\langle\ell_\delta^i\rangle = 0.4253$, $\sigma^i = 0.6458$, and $\langle\ell_\delta^o\rangle = 1.6174$, $\sigma^o = 1.9631$, so that the relative statistical errors are $\frac{\sigma^i}{\langle\ell_\delta^i\rangle\sqrt{N}} \approx 0.15\%$, and $\frac{\sigma^o}{\langle\ell_\delta^o\rangle\sqrt{N}} \approx 0.12\%$.</caption><thead><tr><th>Method</th><th>ε</th><th>CPU time (s)</th><th>⟨ℓ̂<sup>i</sup>⟩<sub>emp</sub></th><th>Rel. Error</th><th>⟨ℓ̂<sup>o</sup>⟩<sub>emp</sub></th><th>Rel. Error</th></tr></thead><tbody><tr><td rowspan="3">SIM1</td><td>10<sup>-1</sup></td><td>1.26 × 10<sup>1</sup></td><td>0.3346</td><td>-21.33%</td><td>1.5072</td><td>-6.81%</td></tr><tr><td>10<sup>-2</sup></td><td>1.39 × 10<sup>2</sup></td><td>0.4385</td><td>3.09%</td><td>1.6888</td><td>4.41%</td></tr><tr><td>10<sup>-3</sup></td><td>1.34 × 10<sup>2</sup></td><td>0.4496</td><td>5.71%</td><td>1.7145</td><td>6.00%</td></tr><tr><td rowspan="3">SIM2</td><td>10<sup>-1</sup></td><td>1.06 × 10<sup>1</sup></td><td>0.4536</td><td>6.65%</td><td>1.4994</td><td>-7.30%</td></tr><tr><td>10<sup>-2</sup></td><td>1.03 × 10<sup>2</sup></td><td>0.4221</td><td>-0.75%</td><td>1.5825</td><td>-2.16%</td></tr><tr><td>10<sup>-3</sup></td><td>1.01 × 10<sup>3</sup></td><td>0.4196</td><td>-1.34%</td><td>1.5961</td><td>-1.32%</td></tr><tr><td rowspan="3">FLA</td><td>10<sup>-1</sup></td><td>1.80 × 10<sup>1</sup></td><td>0.4656</td><td>9.46%</td><td>1.4762</td><td>-8.73%</td></tr><tr><td>10<sup>-2</sup></td><td>1.88 × 10<sup>2</sup></td><td>0.4282</td><td>0.69%</td><td>1.6011</td><td>-1.01%</td></tr><tr><td>10<sup>-3</sup></td><td>1.78 × 10<sup>2</sup></td><td>0.4259</td><td>0.15%</td><td>1.6159</td><td>-0.10%</td></tr><tr><td rowspan="3">CLA</td><td>10<sup>-1</sup></td><td>1.90 × 10<sup>1</sup></td><td>0.4263</td><td>0.24%</td><td>1.6258</td><td>0.52%</td></tr><tr><td>10<sup>-2</sup></td><td>1.86 × 10<sup>2</sup></td><td>0.4253</td><td>-0.0002%</td><td>1.6160</td><td>-0.09%</td></tr><tr><td>10<sup>-3</sup></td><td>1.73 × 10<sup>3</sup></td><td>0.4252</td><td>-0.03%</td><td>1.6174</td><td>-0.004%</td></tr></tbody></table>

(ii) for the boundary local time $\ell_\delta^o$ on the outer sphere:

$$ \pi_0 = 1 - \frac{L}{r_0} \frac{\alpha R \cosh(\alpha(r_0 - R)) + \sinh(\alpha(r_0 - R))}{\alpha R \cosh(\alpha(L - R)) + \sinh(\alpha(L - R))}, \qquad (\text{E.5a}) $$

$$ \mu_0^{(p)} = \frac{(\alpha^2 LR - 1) \sinh(\alpha(L - R)) + \alpha(L - R) \cosh(\alpha(L - R))}{L(\alpha R \cosh(\alpha(L - R)) + \sinh(\alpha(L - R)))}. \qquad (\text{E.5b}) $$

Table E.7 and Figure E.11 present the comparison between three methods, with the same trends as for the circular annulus.

## References

[1] P. Lévy, *Processus Stochastiques et Mouvement Brownien* (Paris, Gauthier-Villard, 1965).

[2] K. Ito and H. P. McKean, *Diffusion Processes and Their Sample Paths* (Springer-Verlag, Berlin, 1965).

[3] M. Freidlin, *Functional Integration and Partial Differential Equations* (Annals of Mathematics Studies, Princeton University Press, Princeton, New Jersey, 1985).

[4] D. S. Grebenkov, *Paradigm Shift in Diffusion-Mediated Surface Phenomena*, Phys. Rev. Lett. **125**, 078102 (2020).

[5] H. P. McKean, *Brownian local time*, Adv. Math. **15**, 91-111 (1975).

[6] A. N. Borodin and P. Salminen, *Handbook of Brownian Motion: Facts and Formulae* (Birkhauser Verlag, Basel-Boston-Berlin, 1996).

[7] S. N. Majumdar, *Brownian functionals in physics and computer science*, Curr. Sci. **89**, 2076 (2007).
<figure><img src="image_12.png" /><figcaption>Figure E.11: Ratio $\Psi_{\text{emp}}(\ell, p|x_0)$ over $\Psi(\ell, p|x_0)$ as a function of $\ell/\langle\ell_\delta\rangle$ for the spherical shell (Fig. 3e), with $p=1$, $R=1$, $L=2$, and $r_0 = 1.5$. The empirical PDF $\Psi_{\text{emp}}$ was estimated from $N = 10^6$ simulated values of $\ell_\delta$ by (a,b) SIM1, (c,d) FLA, (e,f) CLA, with different choices of the boundary layer width: $\varepsilon = 10^{-1}$ (triangles), $\varepsilon = 10^{-2}$ (circles), and $\varepsilon = 10^{-3}$ (squares). (a,c,e) refer to $\ell_\delta^i$ on the inner sphere and (b,d,f) refer to $\ell_\delta^o$ on the outer sphere. The trajectories that never touched the boundary and thus yielded $\ell_\delta = 0$, were excluded.</figcaption></figure>

[8] D. S. Grebenkov, Statistics of boundary encounters by a particle diffusing outside a compact planar domain, J. Phys. A.: Math. Theor. 54, 015003 (2021).

[9] D. S. Grebenkov, An encounter-based approach for restricted diffusion with a gradient drift, J. Phys. A: Math. Theor. 55, 045203 (2022).

[10] P. C. Bressloff, Narrow capture problem: An encounter-based approach to partially reactive targets, Phys. Rev. E 105, 034141 (2022).

[11] D. S. Grebenkov, Statistics of diffusive encounters with a small target: Three complementary approaches, J. Stat. Mech. 083205 (2022).

[12] P. C. Bressloff, Diffusion-mediated surface reactions and stochastic resetting, J. Phys. A: Math. Theor. 55, 275002 (2022).

[13] Z. Benkhadaj and D. S. Grebenkov, Encounter-based approach to diffusion with resetting, Phys. Rev. E 106, 044121 (2022).

[14] D. S. Grebenkov, Depletion of Resources by a Population of Diffusing Species, Phys. Rev. E 105, 054402 (2022).

[15] D. S. Grebenkov, Diffusion-controlled reactions with non-Markovian binding/unbinding kinetics, J. Chem. Phys. 158, 214111 (2023).

[16] D. S. Grebenkov, Adsorption and permeation events in molecular diffusion, Molecules 29, 5012 (2024).

[17] P. C. Bressloff, A probabilistic model of diffusion through a semipermeable barrier, Proc. Roy. Soc. A 478, 20220615 (2022).

[18] P. C. Bressloff, Renewal equation for single-particle diffusion through a semipermeable interface, Phys. Rev. E. 107, 014110 (2023).

[19] P. C. Bressloff, Renewal equations for single-particle diffusion in multilayered media, SIAM J. Appl. Math. 83, 1518-1545 (2023).

[20] E. Fieremans, D. S. Novikov, J. H. Jensen, and J. A. Helpern, Monte Carlo study of a two-compartment exchange model of diffusion, NMR Biomed. 23, 711-724 (2010).

[21] A. Lejay, The snapping out Brownian motion, Ann. Appl. Probab. 26, 1727-1742 (2016).

[22] A. Lejay, Monte Carlo estimation of the mean residence time in cells surrounded by thin layers, Math. Comput. Simul. 143, 65-77 (2018).

[23] D. S. Grebenkov, Probability distribution of the boundary local time of reflected Brownian motion in Euclidean domains, Phys. Rev. E 100, 062110 (2019).

[24] D. S. Grebenkov, Surface Hopping Propagator: An Alternative Approach to Diffusion-Influenced Reactions, Phys. Rev. E 102, 032125 (2020).

[25] K. K. Sabelfeld, *Monte Carlo Methods in Boundary Value Problems* (Springer-Verlag: New York - Heidelberg, Berlin, 1991).

[26] K. K. Sabelfeld and N. A. Simonov, *Random Walks on Boundary for Solving PDEs* (Utrecht, The Netherlands, 1994).

[27] G. N. Milshtein, *Numerical Integration of Stochastic Differential Equations* (Kluwer, Dordrecht, the Netherlands, 1995).

[28] S. Litwin, Monte Carlo simulation of particle adsorption rates at high cell concentration, Biophys. J. 31, 271-277 (1980).

[29] S. Torquato and I. C. Kim, Efficient simulation technique to compute effective properties of heterogeneous media, Appl. Phys. Lett. 55, 1847 (1989).

[30] L. H. Zheng and Y. C. Chiew, Computer simulation of diffusion-controlled reactions in dispersions of spherical sinks, J. Chem. Phys. 90, 322-327 (1989).

[31] S. B. Lee, I. C. Kim, C. A. Miller, and S. Torquato, Random-walk simulation of diffusion-controlled processes among static traps, Phys. Rev. B 39, 11833-11839 (1989).

[32] P. Ossadnik, Multiscaling Analysis of Large-Scale Off-Lattice DLA, Physica A 176, 454 (1991).

[33] L. Batsilas, A. M. Berezhkovskii, and S. Y. Shvartsman, Stochastic model of autocrine and paracrine signals in cell culture assays, Biophys. J. 85, 3659-3665 (2003).
[34] L. Dagdug, A. Berezhkovskii, S. M. Bezrukov, and G. H. Weiss, Diffusion-controlled reactions with a binding site hidden in a channel, J. Chem. Phys. **118**, 2367-2373 (2003).

[35] D. S. Grebenkov, What Makes a Boundary Less Accessible, Phys. Rev. Lett. **95**, 200602 (2005).

[36] D. S. Grebenkov, A. A. Lebedev, M. Filoche and B. Sapoval, Multifractal Properties of the Harmonic Measure on Koch Boundaries in Two and Three Dimensions, Phys. Rev. E **71**, 056121 (2005).

[37] D. S. Grebenkov, Scaling Properties of the Spread Harmonic Measures, Fractals **14**, 231-243 (2006).

[38] P. Levitz, D. S. Grebenkov, M. Zinsmeister, K. Kolwankar and B. Sapoval, Brownian flights over a fractal nest and first passage statistics on irregular surfaces, Phys. Rev. Lett. **96**, 180601 (2006).

[39] M. Deaconu and A. Lejay, A Random Walk on Rectangles Algorithm, Method. Comput. Appl. Probab. **8**, 135-151 (2006).

[40] T. Opplestrup, V. V. Bulatov, G. H. Gilmer, M. H. Kalos, and B. Sadigh, First-passage Monte Carlo algorithm: Diffusion without all the hops, Phys. Rev. Lett. **97**, 230602 (2006).

[41] M. Hall and D. C. Alexander, Convergence and Parameter Choice for Monte-Carlo Simulations of Diffusion MRI, IEEE Trans. Med. Imag. **28**, 1354 (2009).

[42] S. Zein, A. Lejay, and M. Deaconu, An efficient algorithm to simulate a Brownian motion over irregular domains, Commun. Comput. Phys. **8**, 901-916 (2010).

[43] A. M. Berezhkovskii, L. Dagdug, M.-V. Vazquez, V. A. Lizunov, J. Zimmerberg, and S. M. Bezrukov, Trapping of diffusing particles by clusters of absorbing disks on a reflecting wall with disk centers on sites of a square lattice, J. Chem. Phys. **138**, 064105 (2013).

[44] S. K. Ghosh, A. G. Cherstvy, and R. Metzler, Non-universal tracer diffusion in crowded media of non-inert obstacles, Phys. Chem. Chem. Phys. **17**, 1847-1858 (2015).

[45] S. K. Ghosh, A. G. Cherstvy, D. S. Grebenkov, and R. Metzler, Anomalous, non-Gaussian tracer diffusion in heterogeneously crowded environments, New J. Phys. **18**, 013027 (2016).

[46] A. J. Bernoff, A. E. Lindsay, and D. D. Schmidt, Boundary Homogenization and Capture Time Distributions of Semipermeable Membranes with Periodic Patterns of Reactive Sites, Multiscale Model. Simul. **16**, 1411-1447 (2018).

[47] M. Palombo, D. C. Alexander, and H. Zhang, A generative model of realistic brain cells with application to numerical simulation of the diffusion-weighted MR signal, NeuroImage **188**, 391-402 (2019).

[48] A. Ianus, D. C. Alexander, H. Zhang, and M. Palombo, Mapping complex cell morphology in the grey matter with double diffusion encoding MR: A simulation study, NeuroImage **241**, 118424 (2021).

[49] F. Le Vot, S. B. Yuste, E. Abad, and D. S. Grebenkov, First-encounter time of two diffusing particles in two- and three-dimensional confinement, Phys. Rev. E **105**, 044119 (2022).

[50] J. Cherry, A.E. Lindsay, A. Navarro Hernández, and B. Quaife, Trapping of planar Brownian motion: full first passage time distributions by kinetic Monte Carlo, asymptotic, and boundary integral methods, Multiscale Model. Simul. **20**, 1284-1314 (2022).

[51] M. Bossy, E. Gobet, and D. Talay, A symmetrized Euler scheme for an efficient approximation of reflected diffusions, J. Appl. Probab. **41**, 877-889 (2004).

[52] D. S. Grebenkov, Residence times and other functionals of reflected Brownian motion, Phys. Rev. E **76**, 041139 (2007).

[53] M. E. Muller, Some Continuous Monte Carlo Methods for the Dirichlet Problem, Ann. Math. Statist. **27**, 569 (1956).

[54] Y. Zhou, W. Cai, and E. Hsu, Computation of the local time of reflecting Brownian motion and the probabilistic representation of the Neumann problem, Comm. Math. Sci. **15**, 237-259 (2017).

[55] R. D. Schumm and P. C. Bressloff, A numerical method for solving snapping out Brownian motion in 2D bounded domains, J. Comput. Phys. **493**, 112479 (2023).

[56] Y. Scher, S. Reuveni, and D. S. Grebenkov, Escape of a sticky particle, Phys. Rev. Research **5**, 043196 (2023).

[57] Y. Scher, S. Reuveni, and D. S. Grebenkov, Escape from textured adsorbing surfaces, J. Chem. Phys. **160**, 184105 (2024).

[58] D. S. Grebenkov, Encounter-based approach to the escape problem, Phys. Rev. E **107**, 044105 (2023).

[59] I. Binder and M. Braverman, The rate of convergence of the walk of sphere algorithm, Geom. Func. Anal. **22**, 558-587 (2012).

[60] M. Levitin, D. Mangoubi, and I. Polterovich, Topics in Spectral Geometry (Vol. 237), American Mathematical Society (2023).

[61] T. Carlsson, T. Ekholm, and C. Elvingson, Algorithm for generating a Brownian motion on a sphere, J. Phys. A: Math. Theor. **43**, 505001 (2010).

[62] K. Burrage, P. M. Burrage, and G. Lythe, Effective numerical methods for simulating diffusion on a spherical surface in three dimensions, Numer. Algor. **91**, 1577-1596 (2022).

[63] S. B. Yuste, E. Abad, and K. Lindenberg, Exploration and trapping of mortal random walkers, Phys. Rev. Lett. **110**, 220603 (2013).

[64] B. Meerson and S. Redner, Mortality, redundancy, and diversity in stochastic search, Phys. Rev. Lett. **114**, 198101 (2015).

[65] D. S. Grebenkov and J.-F. Rupprecht, The escape problem for mortal walkers, J. Chem. Phys. **146**, 084106 (2017).

[66] S. Maire and E. Tanré, Monte Carlo approximation of the Neumann problem, Monte Carlo Methods Appl. **19**, 201-236 (2013).

[67] D. S. Grebenkov, A fast random walk algorithm for computing the pulsed-gradient spin-echo signal in multiscale porous media, J. Magn. Reson. **208**, 243-255 (2011).

[68] A. Chaigneau and D. S. Grebenkov, A numerical study of the Dirichlet-to-Neumann operator in planar domains, J. Phys. A: Math. Theor. **57**, 445201 (2024).

[69] D. S. Grebenkov, Joint distribution of multiple boundary local times and related first-passage time problems with multiple targets, J. Stat. Mech. 103205 (2020).

[70] C. Costantini, B. Pacchiarotti, and F. Sartoretto, Numerical approximation for functionals of reflecting diffusion processes, SIAM J. Appl. Math. **58**, 73-102 (1998).

[71] G. N. Mil'shtein, Application of the numerical integration of stochastic equations for the solution of boundary-value problems with Neumann's boundary conditions, Theory Probab. Appl. **41**, 170-177 (1996).

[72] G. N. Milstein, Weak approximation of a diffusion process in a bounded domain, Stoch. Stoch. Rep. **62**, 147-200 (1997).

[73] L. Slomiński, Euler's approximations of solutions of SDEs with reflecting boundary, Stoch. Proc. Appl. **94**, 317-337 (2001).

[74] G. N. Milstein and M. V. Tretyakov, *Stochastic Numerics for Mathematical Physics. Scientific Computation* (Springer, Berlin, 2004).
[75] F. Bernal and J. A. Acebrón, A comparison of higher-order weak numerical schemes for stochastic differential equations in bounded domains, Commun. Comput. Phys. 20, 703-732 (2016).

[76] F. Bernal, An implementation of Milstein's method for general bounded diffusions, J. Sci. Comp. 79, 867-890 (2019).

[77] B. Leimkuhler, A. Sharma, and M. V. Tretyakov, Simplest random walk for approximating Robin boundary value problems and ergodic limits of reflected diffusions, Ann. Appl. Probab. 33, 1904-1960 (2023).

[78] D. S. Grebenkov, Analytical representations of the spread harmonic measure density, Phys. Rev. E. 91(5), 052108 (2015).

[79] S. Redner, *A Guide to First Passage Processes* (Cambridge University press, 2001).

[80] Z. Schuss, *Brownian Dynamics at Boundaries and Interfaces in Physics, Chemistry and Biology* (Springer, New York, 2013).

[81] R. Metzler, G. Oshanin, and S. Redner (Eds), *First-Passage Phenomena and Their Applications* (World Scientific Press, Singapore, 2014).

[82] J. Masoliver, *Random Processes: First-passage and Escape* (World Scientific, 2018).

[83] K. Lindenberg, R. Metzler, and G. Oshanin (Eds), *Chemical Kinetics: Beyond the Textbook* (World Scientific, New Jersey, 2019).

[84] L. Dagdug, J. Peña, and I. Pompa-García, *Diffusion Under Confinement: A Journey Through Counterintuition* (Springer, 2024).

[85] D. S. Grebenkov, R. Metzler, and G. Oshanin (Eds), *Target Search Problems* (Springer, Cham, 2024).

[86] E. Hsu, Probabilistic approach to the neumann problem, Commun. Pure Appl. Math. 38, 445 (1985).

[87] V. G. Papanicolaou, The probabilistic solution of the third boundary value problem for second order elliptic equations, Probab. Theory Relat. Fields 87, 27 (1990).

[88] Y. Zhou and W. Cai, Numerical Solution of the Robin Problem of Laplace Equations with a Feynman-Kac Formula and Reflecting Brownian Motions, J. Scient. Comp. 69, 107-121 (2016).

[89] G. A. Brosamler, A probabilistic solution of the Neumann problem, Mathematica Scandinavica, 38, 137-147 (1976).

[90] A. Benchérif-Madani and É. Pardoux, A probabilistic formula for a Poisson equation with Neumann boundary condition, Stoch. Anal. Appl. 27, 739-746 (2009).

[91] P. L. Lions and A. S. Sznitman, Stochastic differential equations with reflecting boundary conditions, Comm. Pure Appl. Math. 37, 511-537 (1984).

[92] J. P. Morillon, Numerical solutions of linear mixed boundary value problems using stochastic representations, Int. J. Numer. Meth. Engng. 40, 387-405 (1997).

[93] D. S. Grebenkov, Efficient Monte Carlo methods for simulating diffusion-reaction processes in complex systems, in "First-Passage Phenomena and Their Applications", Eds. R. Metzler, G. Oshanin, S. Redner (World Scientific Press, Singapore, 2014).

[94] D. S. Grebenkov and S. D. Traytak, Semi-analytical computation of Laplacian Green functions in three-dimensional domains with disconnected spherical boundaries, J. Comput. Phys. 379, 91-117 (2017).