# Boundary local time on wedges and prefractal curves

Yilin Ye, Denis S Grebenkov

► To cite this version:

Yilin Ye, Denis S Grebenkov. Boundary local time on wedges and prefractal curves. Journal of Physics A: Mathematical and Theoretical, 2025, 58 (34), pp.345002. 10.1088/1751-8121/adfb36 hal-05344115

HAL Id: hal-05344115

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

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
# Boundary local time on wedges and prefractal curves

Yilin Ye, Denis S. Grebenkov

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

E-mail: yilin.ye@polytechnique.edu, denis.grebenkov@polytechnique.edu

**Abstract.** We investigate the boundary local time on polygonal boundaries such as finite generations of the Koch snowflake. To reveal the role of angles, we first focus on wedges and obtain the mean boundary local time, its variance, and the asymptotic behavior of its distribution. Moreover, we establish the coupled partial differential equations for higher-order moments. Next, we propose an efficient multi-scale Monte Carlo approach to simulate the boundary local time, as well as the escape duration and position of the associated reaction event on a polygonal boundary. This numerical approach combines the walk-on-spheres algorithm in the bulk with an approximate solution of the escape problem from a sector. We apply it to investigate how the statistics of the boundary local time depends on the geometric complexity of the Koch snowflake. Eventual applications to diffusion-controlled reactions on partially reactive boundaries, including the asymptotic behavior of the survival probability, are discussed.

**Keywords:** boundary local time, diffusion-controlled reactions, first-passage time, survival probability, encounter-based approach, Monte Carlo simulations, fractals
Submitted to:

## 1. Introduction

The boundary local time $l_t$ plays a significant role in the theory of stochastic processes [1–3]. For example, the reflected Brownian motion $X_t$ inside a given Euclidean domain $\Omega \subset \mathbb{R}^d$ with a smooth boundary $\partial\Omega$ satisfies the Skorokhod stochastic equation

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

where $W_t$ is a standard Wiener process in $\mathbb{R}^d$, $D$ is the constant diffusion coefficient, $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 non-decreasing stochastic process that increments only at encounters of $X_t$ with the boundary. The first term in Eq. (1) represents the ordinary Brownian motion
inside $\Omega$, whereas the second term assures its normal reflection back into $\Omega$ after each collision [4–7]. Curiously, the single stochastic equation (1) determines simultaneously two coupled stochastic processes: the position $\mathbf{X}_t$ and the boundary local time $\ell_t$. In addition, the boundary local time can be expressed as

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

where $|\mathbf{x} - \partial\Omega|$ is the Euclidean distance between a point $\mathbf{x}$ and the boundary $\partial\Omega$, and $\Theta(x)$ is the Heaviside step function. Alternatively, one has

$$ \ell_t = \lim_{\epsilon \to 0} \epsilon \mathcal{N}_t^{(\epsilon)}, \quad (3) $$

where $\mathcal{N}_t^{(\epsilon)}$ refers to the number of crossings of a thin boundary layer $\partial\Omega_\epsilon = \{\mathbf{x} \in \Omega : |\mathbf{x} - \partial\Omega| < \epsilon\}$ of width $\epsilon$ near $\partial\Omega$ up to time $t$. In this way, the boundary local time $\ell_t$ can be regarded either as the rescaled residence time in a thin boundary layer or as the rescaled number of encounters with the boundary (despite its name, $\ell_t$ has units of length).

The boundary local time also plays the central role in the encounter-based approach to diffusion-controlled reactions by allowing one to consider general surface reaction or permeation mechanisms [7–17]. There are numerous examples of diffusion-controlled reactions in porous media or on rough catalytic surfaces, including oxygen capture on alveolar surface in the lungs [18–21], heterogeneous catalysis [22–27], electrochemical systems [28–32], and biomolecule recognition [33–35]. As these phenomena are often limited by the diffusive transport of reactants toward a reactive boundary $\partial\Omega$, the geometric complexity profoundly influences the encounter statistics, e.g., via the available contact area in a restricted volume. A key question arises when the boundary exhibits partial reactivity [36–42]: how do surface irregularities modulate the reaction kinetics [42–47]? As a proxy of the time spent near the boundary or the number of collisions, the boundary local time $\ell_t$ on a microscopically rough or fractal-like boundary can help to reveal the relations between diffusion, geometry, and reactivity, linking stochastic trajectories with macroscopic observables.

Despite its theoretical significance, the relationship between the distribution of $\ell_t$ and the geometric form of the boundary remains poorly understood. Although an analytical solution is available for simple geometries like half-spaces (see below), extensions to wedges and polygonal boundaries confront inherent mathematical challenges [48–53]. For instance, an exact solution for the survival probability in a wedge with partially reactive rays is yet unknown. Even on the numerical side, the most advanced algorithms that combine the walk-on-spheres (WOS) method [54] and Skorokhod integral representation [55,56] or the escape-from-a-layer approach [57], allow for simulating the reflected Brownian motion and the boundary local time in domains with smooth boundaries. In turn, undefined or multiple reflection directions near geometric singularities (e.g., corners or cusps) complicate both theoretical formulations
<figure><img src="image_2.png" /><figcaption>Figure 1: (a) Starting from the zeroth generation $g = 0$ (an equilateral triangle with the edges of length $L = 2$), one constructs finite generations $\Omega_g$ of the Koch snowflake iteratively. In the second generation $\Omega_2$, we indicate two starting points used for validation purposes (see Sec. 3.3): at $(0, 5\sqrt{3}/12)$ (red diamond) and at $(1/5, -\sqrt{3}/12)$ (violet square). (b) Local environment near the bottom vertex (red dot) of a finite generation of the Koch snowflake of angle $\alpha = \pi/3$.</figcaption></figure>

and numerical simulations of reflected diffusion in multiscale or self-similar structures like a Koch snowflake (Fig. 1a). This gap hinders the analysis of reaction kinetics on rough boundaries which are often modeled as self-similar fractals, where the statistics of $l_t$ may exhibit geometry-driven anomalies like log-periodic oscillations [58, 59]. In this paper, we develop the escape-from-a-sector (EFS) approach to bridge this gap and to investigate the boundary local time distributions on prefractal curves. We also derive several theoretical results on $l_t$ for wedges that will guide our interpretations of numerical results for the Koch snowflake.

The paper is organized as follows. In Sec. 2, we analyze the boundary local time distribution in wedges, by focusing on the mean value, the variance, and the asymptotic behavior at large $\ell$ (see also Appendix A). Section 3 furnishes the escape-from-a-sector approach for simulating the boundary local time in polygonal domains, and its validation in wedges (Sec. 3.2), and in the Koch snowflake (Sec. 3.3). In Sec. 4, numerical simulations are performed to address different aspects: the mean boundary local time for a wedge with two channels (Sec. 4.1), the boundary local time distribution in the Koch snowflake (Sec. 4.2), and the local persistence exponent (LPE) of the survival probability in the Koch snowflake (Sec. 4.3). Finally, Sec. 5 presents further improvements, eventual applications, and conclusions. Many technical derivations are relegated to Appendices.

## 2. Theoretical results for wedges

In this Section, we present our main theoretical results for a wedge of angle $\alpha$: $\Omega = \{(r, \theta) : r > 0, 0 < \theta < \alpha\}$, written in polar coordinates $(r, \theta)$. Despite the simplicity of this shape, very little is known about the distribution of the boundary local time. As shown in [60], the moment-generating function of $l_t$ is given by

$$S_q(t|x_0) = \langle e^{-q\ell_t} \rangle = \int_0^\infty d\ell \, e^{-q\ell} \rho_\alpha(\ell, t|x_0), \quad q > 0, \quad (4)$$
where $\boldsymbol{x}_0 = (x_0, y_0)$ is the starting point, $\rho_\alpha(\ell, t|\boldsymbol{x}_0)$ is the probability density of $\ell_t$, and $S_q(t|\boldsymbol{x}_0)$ is the survival probability of a particle in the presence of partially reactive boundary with the reactivity parameter $q$. The survival probability satisfies:

$$ 
\begin{aligned}
\partial_t S_q(t|\boldsymbol{x}_0) &= D\Delta S_q(t|\boldsymbol{x}_0) && \text{in } \Omega, \\
\partial_n S_q + qS_q &= 0 && \text{on } \partial\Omega, \\
S_q(0|\boldsymbol{x}_0) &= 1, &&
\end{aligned}
\qquad (5) $$

where $\Delta$ is the Laplace operator, and $\partial_n$ is the normal derivative oriented outwards $\Omega$. In the case $q = \infty$ (i.e., for the Dirichlet boundary condition on $\partial\Omega$), the survival probability $S_\infty(t|\boldsymbol{x}_0)$ can be found explicitly via separation of variables [61–63]. In contrast, the Robin boundary condition does not allow for such a separation in polar coordinates, and we are not aware of such results for partially reactive wedges with $0 < q < +\infty$.

Two exceptions are the upper half-plane (the wedge of angle $\pi$) and the positive quadrant (the wedge of angle $\pi/2$), for which the separation of variables can be realized directly in Cartesian coordinates. In the first case, $\ell_t$ is simply the boundary local time of the reflected Brownian motion on the half-line, whose distribution is well known [5]

$$ \rho_{\pi}(\ell, t|0) = \frac{e^{-\ell^2/(4Dt)}}{\sqrt{\pi Dt}}. \qquad (6) $$

Here and throughout this section, the starting point $\boldsymbol{x}_0$ is located in the origin, i.e., in the vertex of the wedge, in order to highlight the effect of the corner (see Appendix B for a general location of the starting point). In turn, for the positive quadrant, the two rays of the boundary are perpendicular to each other, so that the boundary local time $\ell_t$ is the sum of two independent, identically distributed boundary local times $\ell_t^x$ and $\ell_t^y$ on the horizontal and vertical rays, each obeying the distribution in Eq. (6). As a consequence,

$$ 
\begin{aligned}
\rho_{\pi/2}(\ell, t|0) &= \int_0^\infty d\ell_1 \frac{e^{-\ell_1^2/(4Dt)}}{\sqrt{\pi Dt}} \int_0^\infty d\ell_2 \frac{e^{-\ell_2^2/(4Dt)}}{\sqrt{\pi Dt}} \delta(\ell_1 + \ell_2 - \ell) \\
&= \sqrt{\frac{2}{\pi Dt}} e^{-\ell^2/(8Dt)} \operatorname{erf}\left(\frac{\ell}{\sqrt{8Dt}}\right) && (7)
\end{aligned}
 $$

(see Appendix B for an arbitrary starting point). Note that the corresponding survival probabilities are

$$ S_q(t|0) = \operatorname{erfcx}(q\sqrt{Dt}) \quad (\alpha = \pi), \qquad (8) $$

$$ S_q(t|0) = [\operatorname{erfcx}(q\sqrt{Dt})]^2 \quad (\alpha = \pi/2), \qquad (9) $$

where $\operatorname{erfcx}(z) = e^{z^2}\operatorname{erfc}(z)$ is the scaled complementary error function. For these two cases, it is immediate to compute the moments of the boundary local time; in particular,
one gets the mean and the variance

$$ \langle \ell_t \rangle_0 = \frac{2}{\alpha} \sqrt{\pi Dt}, \qquad (10) $$

$$ \mathrm{Var}_0\{\ell_t\} = \frac{2\pi Dt}{\alpha} \left(1 - \frac{2}{\pi}\right), \qquad (11) $$

where the subscript 0 highlights the starting point at the origin. We aim at extending these results to other wedges.

## 2.1. Mean value

We provide two alternative ways to get the mean boundary local time $\langle \ell_t \rangle_0$: (i) an intuitive symmetry argument and (ii) a direct computation. In the first approach, we consider the starting point $x_0$ to be at the vertex. Let us first assume that $\alpha = 2\pi/n$ with an integer $n$ so that the wedge can be replicated $n$ times to cover the whole plane. Due to normal reflections on the boundary of the wedge, random trajectories of ordinary Brownian motion in the plane are statistically equivalent to those of reflected Brownian motion inside the wedge. The boundary local time $\ell_t$ can thus be represented as the sum of $n$ boundary local times $\ell_t^1, \dots, \ell_t^n$ spent on the rays. As a consequence, the mean value is simply

$$ \langle \ell_t \rangle_0 = n \langle \ell_t^1 \rangle_0 = n \sqrt{\frac{Dt}{\pi}} = \frac{2\pi}{\alpha} \sqrt{\frac{Dt}{\pi}}, \qquad (12) $$

where we used that $\langle \ell_t^1 \rangle_0 = \sqrt{Dt/\pi}$. Note that this mean is twice smaller than the mean boundary local time from Eq. (10) on the whole horizontal line, which can be seen as being composed of two rays. We emphasize that this simple argument does not allow one to get higher-order moments, nor the distribution of $\ell_t$, because $\ell_t^1, \dots, \ell_t^n$ are not independent. Indeed, if planar Brownian motion spends more time near one ray, it generally spends less time near the other rays.

In the next step, one can consider $\alpha = 2\pi m/n$ with integer $m$ and $n$. Replicating this wedge $n$ times, one covers the whole plane $m$ times, as if the plane was replicated $m$ times. One can think that Brownian motion switches between $m$ parallel planar “layers”, thus spending on average $\langle \ell_t^1 \rangle_0/m$ on one ray of one layer. As a consequence, we get again $\langle \ell_t \rangle_0 = n \langle \ell_t^1 \rangle_0/m = (2\pi/\alpha)\sqrt{Dt/\pi}$. Finally, as “rational” angles of the form $2\pi m/n$ are dense in the continuum set of all angles, we conclude that

$$ \langle \ell_t \rangle_0 = \frac{2\pi}{\alpha} \sqrt{\frac{Dt}{\pi}} = \frac{2\sqrt{\pi Dt}}{\alpha} \qquad (13) $$

for any wedge.

The second method relies on the following representation of the mean boundary local time (see [60] and Appendix A)

$$ \langle \ell_t \rangle_{x_0} = \int_0^t dt' \int_{\partial\Omega} d\boldsymbol{x} \, DG_0(\boldsymbol{x}, t'|\boldsymbol{x}_0), \qquad (14) $$
where $G_0(x, t|x_0)$ is the propagator with Neumann boundary condition. In Appendix C, we provide the exact expression for this propagator for any wedge, as well as the evaluation of the double integral in Eq. (14). This “brute-force” computation yields again Eq. (13). Moreover, we also obtain $\langle \ell_t \rangle_{x_0}$ for any starting point $x_0$.

## 2.2. Variance and higher-order moments

To access higher-order moments, we establish in Appendix A a system of coupled partial differential equations between the moments of the boundary local time $\ell_t$. Their solution requires the knowledge of the propagator $G_0(x, t|x_0)$ that we have already used in Eq. (14) for evaluating the mean boundary local time. Even though this propagator is known explicitly for any wedge (see Appendix C.1), the resulting expressions are rather cumbersome even for the second moment of $\ell_t$ (see Eq. (A.10)). For this reason, we do not further explore this general direction and restrict our attention to the starting point $x_0$ located at the origin. In this case, the only length scale of the problem is $\sqrt{Dt}$, and a basic dimensional argument implies that $\langle \ell_t^k \rangle_0 = \Lambda_k(\alpha)(Dt)^{k/2}$, where $\Lambda_k(\alpha)$ is (yet unknown) dimensionless prefactor that depends on $k$ and $\alpha$. In particular, Eq. (11) for the variance, which was valid for $\alpha = \pi$ and $\alpha = \pi/2$, can be generalized as

$$ \mathrm{Var}_0\{\ell_t\} = v_\alpha \frac{2\pi Dt}{\alpha} \left(1 - \frac{2}{\pi}\right), \qquad (15) $$

with a prefactor $v_\alpha$ that depends only on the angle. In Appendix C.3, we compute this prefactor exactly in the case when $\alpha = \pi/n$, with an integer $n$. In particular, we retrieve $v_\pi = v_{\pi/2} = 1$ but show that $v_\alpha$ increases as $\alpha$ decreases.

## 2.3. Asymptotic behavior at large $\ell$

We now inspect the asymptotic behavior of the probability density $\rho_\alpha(\ell, t|x_0)$ of $\ell_t$. For a bounded domain, the Laplace transform of this probability density admits the following spectral expansion [60]:

$$ \tilde{\rho}(\ell, p|x_0) = \int_0^\infty dt \, e^{-pt} \rho(\ell, t|x_0) = \frac{1}{p} \sum_k \mu_k^{(p)} e^{-\mu_k^{(p)} \ell} V_k^{(p)}(x_0) \int_{\partial\Omega} d\mathbf{x} \, V_k^{(p)}(\mathbf{x}), \qquad (16) $$

where $\mu_k^{(p)}$ and $V_k^{(p)}$ are the eigenvalues and eigenfunctions of the generalized Steklov problem:

$$ \begin{aligned} (p - D\Delta)V_k^{(p)} &= 0 && \text{in } \Omega, \\ \partial_n V_k^{(p)} &= \mu_k^{(p)} V_k^{(p)} && \text{on } \partial\Omega. \end{aligned} \qquad (17) $$

This spectral problem is known to have a discrete spectrum [64], so that the eigenmodes can be enumerated by a positive integer $k$. In turn, as the boundary of a wedge is unbounded, the spectrum of the Steklov problem is not discrete anymore. Nevertheless,
for a wedge of angle $\alpha < \pi$, there exists at least one eigenvalue below the essential spectrum (see, e.g., [65–67] and references therein), and this eigenvalue determines the asymptotic behavior at large $p$ or, equivalently, at short times.

To access this behavior, let us inspect the spectrum of the Robin Laplacian in the wedge of angle $\alpha = 2\Phi$: $\Omega' = \{(r, \theta) : |\theta| < \Phi\}$ that we rotated by angle $\Phi$ for convenience:

$$ \begin{aligned} -\Delta u &= \lambda u && \text{in } \Omega', \\ \partial_n u &= \mu u && \text{on } \partial\Omega', \end{aligned} \quad (18) $$

with a prescribed parameter $\mu > 0$. Even though the spectrum is continuous, there may exist a finite number of negative eigenvalues. In particular, the smallest eigenvalue is $\lambda_1 = -\mu^2/\sin^2(\Phi)$, whereas the associated eigenfunction is $u_1 = \exp(-\mu x_0/\sin(\Phi))$ [68]. Setting $\lambda_1 = -p/D$ and employing the duality between the Steklov and Robin problems [64], one sees that there exists a Steklov eigenpair:

$$ \mu_1^{(p)} = \sin(\Phi) \sqrt{p/D}, \quad V_1^{(p)}(x_0, y_0) = C_1 e^{-\mu_1^{(p)} x_0 / \sin(\Phi)}, \quad (19) $$

where $x_0 = (x_0, y_0)$, and $C_1$ is the normalization constant that ensures $L^2(\partial\Omega)$-normalization of $V_1^{(p)}$:

$$ \int_{\partial\Omega'} d\boldsymbol{x} [V_1^{(p)}(\boldsymbol{x})]^2 = 1 \quad \Rightarrow \quad C_1^2 = \frac{\mu_1^{(p)}}{\tan(\Phi)}. \quad (20) $$

As $\mu_1^{(p)}$ is the smallest eigenvalue, it controls the asymptotic behavior of Eq. (16) at large $p$. As a consequence, we have as $p \to \infty$

$$ \tilde{\rho}_{\alpha}(\ell, p|\boldsymbol{x}_0) \approx \frac{1}{p} \mu_1^{(p)} e^{-\mu_1^{(p)} \ell} V_1^{(p)}(\boldsymbol{x}_0) \int_{\partial\Omega'} d\boldsymbol{x} V_1^{(p)}(\boldsymbol{x}) = \frac{2 \sin(\Phi)}{\sqrt{pD}} e^{-\sqrt{p/D}(x_0 + \ell \sin(\Phi))}, \quad (21) $$

from which the inverse Laplace transform yields

$$ \rho_{\alpha}(\ell, t | \boldsymbol{x}_0) \approx \frac{2 \sin(\alpha/2)}{\sqrt{\pi Dt}} e^{-(r_0 \cos(\theta_0 - \alpha/2) + \ell \sin(\alpha/2))^2/(4Dt)}, \quad (22) $$

where we wrote $x_0$ in terms of the polar coordinates in our conventional wedge $\Omega = \{(r, \theta) : 0 < \theta < \alpha\}$. This expression determines the short-time or, equivalently, the large-$\ell$ behavior.

In contrast, there is no isolated eigenvalue for a wedge of angle $\alpha > \pi$, and the bottom of the essential spectrum is $\sqrt{p/D}$. As a consequence, one may expect the asymptotic behavior $\rho_{\alpha}(\ell, t|0) \propto e^{-\ell^2/(4Dt)}$ for any $\alpha > \pi$, with eventual power-law corrections.
## 2.4. Asymptotic behavior at small $\ell$

To complete this section, we briefly discuss the asymptotic behavior of the probability density $\rho_{\alpha}(\ell, t|0)$ at small $\ell$. From a dimensional argument, this is a function of $\ell/\sqrt{Dt}$. According to the explicit solutions in Eqs. (6, 7), one can expect a power-law behavior,

$$ \rho_{\alpha}(\ell, t|0) \propto \frac{1}{\sqrt{Dt}} (\ell/\sqrt{Dt})^{\beta-1} \quad (\ell \to 0), \qquad (23) $$

with an exponent $\beta$ that depends only on the angle $\alpha$. In particular, we have $\beta = 1$ for $\alpha = \pi$ and $\beta = 2$ for $\alpha = \pi/2$. According to Eq. (4), the asymptotic behavior in Eq. (23) is tightly related to the large-$q$ asymptotic behavior of the survival probability $S_q(t|0)$ as

$$ S_q(t|0) \propto (q\sqrt{Dt})^{-\beta} \quad (q \to \infty). \qquad (24) $$

As the survival probability $S_q(t|0)$ is a function of $q\sqrt{Dt}$, its large-$q$ asymptotic behavior is equivalent to the long-time asymptotic behavior. Since we are not aware of earlier theoretical studies on this quantity for partially reactive wedges, we propose the following heuristic argument. At large $q$ and $t$, the asymptotic behavior of the survival probability $S_q(t|0)$ is expected to be close to that of $S_{\infty}(t|r_0)$. Indeed, a particle that managed to survive for a long time should rapidly move away from the wedge boundary and keep avoiding it. The long-time asymptotic behavior of $S_{\infty}(t|r_0)$ is well known (see, e.g., [69–71]): $S_{\infty}(t|r_0) \propto (r_0/\sqrt{Dt})^{\pi/\alpha}$, so that we conjecture that $\beta = \pi/\alpha$. This conjecture agrees with the aforementioned exact values of $\beta$ for $\alpha = \pi$ and $\alpha = \pi/2$. Further numerical analysis of this conjecture will be reported elsewhere.

## 3. Escape-from-a-sector approach

How can one simulate the boundary local time? The simplest method is to discretize the space, such that the reflected Brownian motion is modeled by a random walk on a lattice with a small enough spacing $\epsilon$, which would require plenty of calculation resources for long trajectories. The continuous space simulations can also be realized via the standard walk-on-sphere (WOS) algorithm by Muller [54], combined with constant displacements inside a thin boundary layer [55]. Moreover, Schumm and Bressloff [56] implemented the Skorokhod integral representation in planar bounded domains with smooth boundaries.

However, as the boundary of a polygonal domain is not smooth near vertices, we are not aware of Monte Carlo techniques for simulating efficiently the boundary local time in such settings. To bridge this gap, we develop here the escape-from-a-sector approach based on the spectral decomposition of a suitable escape problem.

We first recall briefly the WOS algorithm for simulating Brownian motion in a confining domain $\Omega \subset \mathbb{R}^2$ with the boundary $\partial\Omega$ [54]. From a given starting point $\boldsymbol{x}_0 \in \Omega$, one draws a disk $B_{\rho}(\boldsymbol{x}_0)$ of radius $\rho = |\boldsymbol{x}_0 - \partial\Omega|$ centered at $\boldsymbol{x}_0$. A continuous trajectory of Brownian motion crosses the boundary of the disk at some random point $\boldsymbol{x}_1$, which is uniformly distributed on $\partial B_{\rho}(\boldsymbol{x}_0)$, at some random escape time $\tau_1$, whose
<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic of a random trajectory of a particle moving in a wedge (with black dots representing sampled random positions). An escape event is initiated when the particle has entered a circular sector of small radius $\epsilon$ and angle $\alpha$. The blue dot indicates the current position $\mathbf{x}_0$ of the particle inside the sector, which in local polar coordinates reads as $(r_0, \theta_0)$. The random escape position $\mathbf{X}_{\tau} = (\epsilon, \theta_{\tau})$ (red dot) is located on the arc $\partial\Omega_D$. The boundary of the circular sector is composed of two reflecting segments $\Gamma$ and the absorbing arc $\partial\Omega_D$.</figcaption></figure>

distribution is known explicitly [70]. As a consequence, a detailed simulation of the Brownian trajectory $\mathbf{X}_t$ inside the disk can be replaced by generating a random escape point $\mathbf{X}_{\tau_1} = \mathbf{x}_1$ at the escape time $\tau_1$. Repeating this procedure, one samples a sequence of points $\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_k$ of a random trajectory at random times $t_k = t_{k-1} + \tau_k$, where $\tau_j$ are independent escape times. This procedure is iterated until the distance to the boundary becomes smaller than a prescribed threshold $\epsilon$, e.g., the width of a boundary layer. If the boundary is perfectly absorbing, the simulation is stopped, and the current time and position are recorded as the first-passage time to the boundary and its location. In our setting, however, the boundary $\partial\Omega$ is not perfectly absorbing, so that one needs to simulate reflections on the boundary after the first arrival into the boundary layer. This is the most time-consuming step of the former approaches.

To overcome this limitation, one can simulate the escape from that layer by generating the escape time $\tau$, the escape position $\mathbf{x}'$, and the boundary local time $\ell_{\tau}$ acquired up to the escape event [57]. These random variables can be generated from their distributions that are known for a flat layer, which can be considered as a local approximation of a smooth boundary. This approach provides an accurate framework for simulating the boundary local time in Euclidean domains with smooth boundaries. However, it is not applicable near corners. In this section, we propose an alternative solution, which allows us to handle polygonal boundaries.

Let us consider such a situation when the current position of the particle is closer to a corner of angle $\alpha$ than a prescribed threshold $\epsilon$ (Fig. 2). After diffusing inside the sector of angle $\alpha$ and radius $\epsilon$, the particle escapes this sector through a random point on the arc $\partial\Omega_D$. An accurate simulation of the Brownian trajectory inside the sector would involve multiple reflections on the segments and could result in accumulated errors. We aim at replacing this time-consuming step by generating a single escape
event (we therefore use the name “escape-from-a-sector” approach). For this purpose, one needs to generate three random variables: the escape time $\tau$, the escape position on the arc (characterized by the angle $\theta_{\tau}$), and the boundary local time $\ell_{\tau}$ acquired up to $\tau$. Even though their joint distribution can be formally found [72], its cumbersome form is not suitable for efficient simulations. For this reason, we propose to substitute random realizations of $\tau$, $\theta_{\tau}$, and $\ell_{\tau}$ by their mean values. In other words, we update the counters upon the escape event as:

$$t_{n+1} = t_n + \mathbb{E}_{x_0}\{\tau\}, \quad \ell_{n+1} = \ell_n + \mathbb{E}_{x_0}\{\ell_{\tau}\}, \quad \mathbf{x}_{n+1} = (\varepsilon, \mathbb{E}_{x_0}\{\theta_{\tau}\}), \quad (25)$$

where the escape position $\mathbf{x}_{n+1}$ is given in local polar coordinates of the sector. We can use the following expressions for the mean values that we derive in Appendix D:

$$\mathbb{E}_{x_0}\{\tau\} = \frac{\varepsilon^2 - r_0^2}{4D}, \quad (26)$$

$$\mathbb{E}_{x_0}\{\ell_{\tau}\} = \frac{2\varepsilon}{\alpha} - r_0 \frac{\cos(\theta_0) + \cos(\alpha - \theta_0)}{\sin(\alpha)} - \frac{2\varepsilon}{\alpha} \sum_{n=1}^{\infty} \frac{(1 + (-1)^n)}{\nu_n^2 - 1} \cos(\nu_n \theta_0) (r_0/\varepsilon)^{\nu_n}, \quad (27)$$

and

$$\mathbb{E}_{x_0}\{\theta_{\tau}\} = \frac{\alpha}{2} \left( 1 - 4 \sum_{n=1}^{\infty} \frac{1 - (-1)^n}{\pi^2 n^2} \cos(\nu_n \theta_0) (r_0/\varepsilon)^{\nu_n} \right), \quad (28)$$

where $\nu_n = \pi n / \alpha$. Since $r_0 < \varepsilon$, both sums converge rapidly.

When one needs to generate the boundary local time $\ell_t$, the simulation is stopped when the time counter $t_n$ exceeds a prescribed time $t$. Repeating simulations $N$ times, one gets the empirical statistics of $\ell_t$, from which its moments and the empirical density (rescaled histogram) can be estimated. In turn, the computation of the survival probability $S_q(t|x_0)$ requires a different stopping condition: a simulated trajectory is stopped when the boundary local time $\ell_t$ exceeds a random threshold $\hat{\ell}$ with the exponential distribution: $\mathbb{P}\{\hat{\ell} > \ell\} = e^{-q\hat{\ell}}$ [7]. In other words, the simulation is stopped when $\ell_n > \hat{\ell}$, and the first-reaction time $T$ is assigned to be $t_n$. Repeating simulations $N$ times, one gets the empirical statistics of $T$, from which its moments, the empirical density, and the survival probability can be estimated.

### 3.1. Implementation

We launch Monte Carlo simulations of the reflected Brownian motion starting from a fixed point $x_0$ inside a given domain, e.g., a wedge or a polygon like finite generations of the Koch snowflake. The latter can be constructed iteratively starting from an equilateral triangle of length $L = 2$ and replacing each segment by a simple generator (Fig. 1a). In this way, the $g$-th generation $\Omega_g$ of the Koch snowflake is a polygonal domain whose boundary is composed of $3 \times 4^g$ segments of length $h_g = L/3^g$. As $g$ increases, one adds finer and finer geometrical details such that the perimeter of the boundary, $3 \times 4^g \cdot h_g$, diverges. The limiting domain $\Omega_{\infty}$ has a fractal boundary that is characterized by the fractal dimension $d_f = \ln 4 / \ln 3 \approx 1.26$ [73].
<figure><img src="image_4.png" /><figcaption>Figure 3: Comparison of empirical PDFs of the boundary local time obtained by Monte Carlo simulations in the Koch snowflake at $g = 10$ with $L = 2$, $D = 1$, $t = 10^{-7}$, and $N = 10^6$ particles. The boundary layer thickness is chosen as $\varepsilon = 10^{-5}$ (blue squares), $\varepsilon = 10^{-6}$ (orange diamonds), and $\varepsilon = 10^{-7}$ (green dots).</figcaption></figure>

A prescribed parameter $\varepsilon$ is introduced as the thickness of the boundary layer. When the current position $\boldsymbol{x}_n$ of the particle is far away from the boundary, i.e., the distance $\rho = |\boldsymbol{x}_n - \partial\Omega|$ is larger than $\varepsilon$, the WOS algorithm is executed. When $\boldsymbol{x}_n$ is near the boundary but far away from corners, i.e., $\rho < \varepsilon$ but $\rho_v = |\boldsymbol{x} - \partial\Omega_v| > \varepsilon$ where $\partial\Omega_v$ is the ensemble of vertices, the EFL method [57] with the flat-layer approximation is exploited. If $\boldsymbol{x}_n$ is near any vertex, i.e., $\rho_v < \varepsilon$, the EFS approach is employed.

The choice of $\varepsilon$ is a compromise between accuracy and rapidity. For an accurate computation, one needs to ensure that $\varepsilon$ is smaller than the smallest geometric feature of the boundary. For instance, for the $g$-th generation of the Koch snowflake, one needs $\varepsilon \ll h_g$. Figure 3 compares three empirical probability density functions (PDFs) of the boundary local time $\ell_t$ in $\Omega_{10}$, obtained with different values of $\varepsilon$. As $h_{10} \approx 3.4 \times 10^{-5}$, both choices $\varepsilon = 10^{-6}$ and $\varepsilon = 10^{-7}$ are small enough and thus yield almost indistinguishable PDFs. In turn, the larger value $\varepsilon = 10^{-5}$ yields a different (wrong) PDF. According to this brief verification, we use $\varepsilon = 10^{-6}$ for generations $g$ up to 10. In turn, $\varepsilon = 10^{-7}$ is taken for higher generation $g = 11, 12$. A home-built code was written in Fortran 90 for numerical simulations, while data were analyzed in MATLAB and Python.

It was checked that the CPU time for a numerical simulation of $\ell_t$ depends on four parameters: the generation $g$, the number of particles $N$, the time $t$, and the boundary layer thickness $\varepsilon$, such that CPU $\propto gNt/\varepsilon$. The dependence of CPU on $g$ is almost linear because we exploit the geometry-adapted fast random walk (GAFRW) algorithm for distance calculations [74]. When the stopping condition concerns the reactivity $q$, we have CPU $\propto gN/(q\varepsilon)$. For small reactivity $q$ (almost inert surfaces), a larger number
<figure><img src="image_5.png" /><figcaption>Figure 4: The coefficient $A_\alpha$ of the large-$\ell$ asymptotic behavior of the probability density distribution $\rho_\alpha(\ell, t|0)$ as a function of angle $\alpha$. Circles present the numerical fitting via Eq. (29) from the empirical density obtained by Monte Carlo simulations with $N = 10^6$ particles and $Dt = 1$ for each angle, whereas the solid line shows the theoretical prediction (30). The width of the boundary layer is taken as $\epsilon = 10^{-3}$.</figcaption></figure>

of particles should be taken to diminish statistical fluctuations.

## 3.2. Validation of EFS approach on wedges

To validate the EFS approach on wedges, we first investigate the statistics of the boundary local time for different wedge angles. For this purpose, we set the starting point at the vertex of the wedge of angle $\alpha$. From simulations, we get an excellent agreement with Eq. (7) for $\alpha = \pi/2$ (figure is not shown). We also retrieve numerically the Gaussian right tail of the distribution:

$$ \rho_{\alpha}(\ell, t|0) \propto e^{-\ell^2/(A_{\alpha}Dt)} \quad (\ell \to \infty), \qquad (29) $$

$$ A_{\alpha} = 4/\sin^2(\alpha/2), \qquad (30) $$

for a broad range of angles from $\pi/20$ to $\pi$ (Fig. 4), in perfect agreement with Eq. (22).

Another validation of the EFS approach employs our theoretical prediction (10) for a broad range of angles from $\pi/20$ to $2\pi$ (Fig. 5a). The relative error between theoretical and numerical values of $\langle\ell_t\rangle_0$ is small and can be attributed to a finite statistics and eventual (minor) errors of our approximate simulation scheme. Figure 5b shows the prefactor $v_\alpha$ in Eq. (15) as a function of the angle $\alpha$ of the wedge. For angles of the form $\alpha = \pi/n$, theoretical values obtained in Appendix C.3 are in good agreement with numerical predictions. Curiously, the factor $v_\alpha$ exhibits a non-monotonous dependence on $\alpha$ and increases as $\alpha \to 0$.
<figure><img src="image_6.png" /></figure>

Figure 5: (a) Relative error $|\langle l_t \rangle_0^{\text{num}} / \langle l_t \rangle_0 - 1|$ of the mean boundary local time $\langle l_t \rangle$ as a function of the angle $\alpha$ of the wedge, with $Dt = 1$. The analytical values are given by Eq. (10), whereas numerical values were obtained by Monte Carlo simulations with $N = 10^6$ particles and $\varepsilon = 10^{-3}$. Positive values are shown as squares ($\alpha = \pi/4$), while negative values are shown as diamonds. (b) The prefactor $v_\alpha$ from Eq. (15) as a function of the angle $\alpha$ of the wedge. The analytical values (red squares) are given by Eq. (C.34). Numerical values (black crosses) were obtained by Monte Carlo simulations with $N = 10^6$ particles and $\varepsilon = 10^{-3}$. The thin dashed line presents 1 (the exact value of $v_\alpha$ for $\alpha = \pi$ and $\alpha = \pi/2$).

### 3.3. Validation via spread harmonic measures

For a partially reactive surface, the spread harmonic measure characterizes the spatial distribution of successful reaction events [75, 76]. In order to validate the accuracy of the proposed EFS method in polygonal domains, we compute numerically the spread harmonic measure for several generations of the Koch snowflake. More precisely, we compute the probability $p_k$ of reacting on the $k$-th segment of the boundary of $\Omega_g$. If the starting point is located at the center of the equilateral triangle $\Omega_0$ or the first generation of the Koch snowflake $\Omega_1$, the symmetry implies equal probabilities $p_k$ for all segments, regardless of the reactivity $q$. We checked this statement numerically (results are not shown). Besides, we examine starting points out of the center of the second generation $\Omega_2$, namely, $(0, 5\sqrt{3}/12)$ and $(1/5, -\sqrt{3}/12)$, as shown in Fig. 1b. As a benchmark, an implementation of a finite-element method (FEM) was realized, where the spectral expansion of the spread harmonic measure was truncated (see [7]) to 30 terms, and the maximal mesh size was 0.005 [77]. The accuracy of this FEM computation was checked by varying the truncation order and the mesh size. In Monte Carlo simulations, we set $\varepsilon = 10^{-3} \ll h_2 \approx 2.2 \times 10^{-1}$ to ensure an accurate estimate of the boundary local time.

Figure 6 shows the probabilities $\{p_k\}$ obtained by FEM (blue lines) and Monte
<figure><img src="image_7.png" /></figure>

Figure 6: Spread harmonic measure on the second generation $\Omega_2$ of the Koch snowflake: comparison between the finite-element method (blue lines) and our Monte Carlo technique (orange dots). Parameters are: the length $L = 2$, the boundary layer width $\varepsilon = 10^{-3}$, particle number $N = 10^6$, reactivity $q \in [10, 1, 0.1]$ for each row from the top to the bottom, and the starting point $x_0$ is located at $(0, 5\sqrt{3}/12)$ for the left column and at $(1/5, -\sqrt{3}/12)$ for the right column. The segment index $k$ starts from the left bottom and increases anti-clockwise.

Carlo (MC) techniques (orange dots). Since the first point $(0, 5\sqrt{3}/12)$ is located on the vertical axis of symmetry, we observe the symmetrical pattern of $\{p_k\}$ in three panels in the left column. For the second point $(1/5, -\sqrt{3}/12)$, the distribution exhibits six peaks as expected. With the decrease of reactivity, from top to bottom, we observe fewer differences in $\{p_k\}$ among segments. Both methods provide almost identical results that demonstrate the accuracy and capacity of the EFS approach.

## 4. Numerical results

In this Section, we present our main numerical results for finite generations $\Omega_g$ of the Koch snowflake, when the starting point is located at one of the vertices. As the zeroth generation $\Omega_0$ is just an equilateral triangle of side $L$, the distribution of the boundary local time in $\Omega_0$ is expected to be close to $\rho_{\pi/3}(\ell, t|0)$ in the wedge of angle $\pi/3$, when $t$ is much smaller than $L^2/D$. Indeed, when $t \ll L^2/D$, only a few trajectories of the
<figure><img src="image_8.png" /><figcaption>Figure 7: Schematic of a wedge of angle $\alpha = \pi/2$ with two additional rectangular channels on each ray. The starting point $\mathbf{x}_0$ is located at the vertex of the wedge (red dot).</figcaption></figure>

reflected Brownian motion can reach the opposite side of the triangle and thus "feel" the difference between the triangle and the wedge. At the next iteration, the edges of the equilateral triangle are "decorated" by small triangles. Is this modification of the boundary "beneficial", i.e., does it enhance or diminish the average number of encounters with the boundary? How does it reshape the distribution of $\ell_t$? Answering these basic questions is actually not simple. On the one hand, the perimeter of the boundary, which is typically accessible to the particle up to time $t$, is increased, that may increase $\langle \ell_t \rangle$. On the other hand, two straight segments of the triangle, which were easily accessible to the diffusing particle, are now removed, while the added longer segments are less accessible due to diffusion screening (see [18, 22, 44, 78–80] and references therein). The iterative construction of the next generations makes the answers even more difficult. For this reason, we first consider a simple setting of a wedge with two channels (Sec. 4.1) and then proceed to finite generations of the Koch snowflake (Sec. 4.2).

### 4.1. Wedge with two channels

Before proceeding to the analysis of the boundary local time in self-similar polygonal domains, we consider a minor alteration of the wedge by adding two identical rectangular channels on both rays (Fig. 7). As in the case of the first generation of the Koch snowflake, it is not clear a priori whether such a modification would increase or decrease the mean boundary local time. We investigate here $\langle \ell_t \rangle_0$ as a function of the parameters $a, b$, and $c$ of the channels (here we restrict our discussion to the case $a = b$).

Figure 8 shows the mean boundary local time in the domain with $c = 5$ and $a = b$ at $t = 10, 100, 400$ (with $D = 1$). At $t = 10$ (shown by triangles), the diffusion length $\sqrt{Dt} \approx 3.3$ is smaller than $c$, i.e., most trajectories do not reach the channel, and $\langle \ell_t \rangle_0$ is very close to that of the wedge of angle $\pi/2$, independently of the channel width $a$. At $t = 100$, the diffusion length $\sqrt{Dt} = 10$ is comparable to $c$ so that the channels start to affect the statistics of encounters; in particular, the mean boundary local time increases and then decreases, with the maximum around $a = 3$. Moreover, when $a \ge 9$, $\langle \ell_t \rangle_0$ becomes smaller than the mean value for the wedge without channels. In other words, we see that the presence of the channels can either increase or decrease the mean
<figure><img src="image_9.png" /><figcaption>Figure 8: Mean boundary local time $\langle l_t \rangle_0$ in a wedge of angle $\pi/2$ with two channels (Fig. 7) as a function of $a=b$, with $c=5$, $D=1$, $N=10^6$ particles and $\varepsilon=10^{-2}$ for different times: $t=10$ (blue triangle), $t=100$ (orange square), and $t=400$ (green dot). Horizontal lines indicated $4\sqrt{Dt}/\pi$ that corresponds to $\langle l_t \rangle_0$ in the wedge without channels (i.e., $a=b=0$).</figcaption></figure>

boundary local time, depending on its size: it is beneficial for $a \lesssim \sqrt{Dt}$ and detrimental for $a \gtrsim \sqrt{Dt}$. This is not surprising. In the former case, the larger accessible perimeter tends to increase $\langle l_t \rangle_0$. In turn, when $a \gtrsim \sqrt{Dt}$, the particle does access the whole channel, and the accessible perimeter is smaller. At $t=400$, the above observation remains valid, but the maximum is shifted to $a \approx 6$. This observation illustrates that the diffusion length scale $\sqrt{Dt}$ is not sufficient to determine the maximum, which depends on the shape of the whole domain.

## 4.2. Statistics of the boundary local time in the Koch snowflake

We employ the EFS method to investigate systematically the statistics of the boundary local time in the Koch snowflake, for different times $t$ and generations $g$. In this study, we aim at revealing the effect of boundary complexity onto the boundary local time. For this reason, we start all simulations from a vertex of angle $\pi/3$ (see Fig. 1b) and consider a broad range of intermediate times, such that $t \ll L^2/D = 4$. In this case, the diffusing particle does not have time to reach the central zone of the Koch snowflake and thus to “feel” confinement in a bounded domain. In other words, the particle diffuses near a prefractal boundary, as if the domain was unbounded. In this way, we can reveal the effect of boundary irregularity as compared to the flat boundary of the wedge of angle $\pi/3$.

At very short times ($t \ll h_g^2/D$), the local environment of the Koch snowflake is identical to the vicinity of the vertex of a wedge of angle $\pi/3$, and the distribution of
<figure><img src="image_10.png" /><figcaption>Figure 9: Comparison of empirical PDFs of the boundary local time $\ell_t$ obtained by Monte Carlo simulations in the Koch snowflake with $L = 2$, $D = 1$, $t = 10^{-7}$, and $N = 10^8$ particles at different generations: $g = 6, 7$ (dots), $g = 8$ (up triangles), $g = 9$ (down triangles), $g = 10$ (diamonds), $g = 11$ (squares), $g = 12$ (circles).</figcaption></figure>

the boundary local time should be identical to that in the wedge. In this regime, the distribution does not depend on $g$, given that the condition $t \ll h_g^2/D$ is satisfied. As soon as $t \gtrsim h_g^2/D$, the particle starts to “feel” irregularities of the local environment, and the distribution of $\ell_t$ starts to depend on $g$. This is confirmed by Fig. 9 that presents empirical PDFs of $\ell_t$ at $t = 10^{-7}$ (and $D = 1$). For generations $g \le 7$, $h_g^2/D = 4 \times 9^{-g}$ remains much larger than $t = 10^{-7}$, and the PDFs are identical to those of the wedge of angle $\pi/3$. In turn, for larger $g$, we observe the effect of the generation on the PDFs: higher generations lead to broader distributions.

As the boundary is self-similar, an appropriate rescaling of time $t$ should keep the distribution of the boundary local time invariant. In fact, let us consider the distribution of $\ell_t$ for a given generation $\Omega_g$. The next generation $\Omega_{g+1}$ differs from $\Omega_g$ by adding the geometric details at the new smallest scale $h_{g+1} = h_g/3$. If time $t$ is reduced by a factor of 9, the particle diffusing in $\Omega_{g+1}$ would effectively explore the environment of $\Omega_g$. In other words, the PDF of $\ell_t/9$ in $\Omega_{g+1}$ is expected to be identical to the PDF of $\ell_t/3$ in $\Omega_g$. This self-similarity is indeed observed in Fig. 10. Note that $\ell_t$ was rescaled by its mean value $\langle \ell_t \rangle_0$.

Let us now inspect the mean value of the boundary local time. A systematic study was done for the Koch snowflake at different times $t$ and generations $g$. In order to compare these cases, we introduce the time scale $t_g = h_g^2/D = L^2/(9^g D)$, and consider $\langle \ell_t \rangle_0$ as a function of $t/t_g$ (Fig. 11). For small $g$ and thus small $t/t_g$, the particle explores the local environment of the wedge, yielding the rescaled mean value $\langle \ell_t \rangle_0/\sqrt{Dt} = 6/\sqrt{\pi} \approx 3.4$, which is shown by a horizontal line. When $t$ becomes comparable to $t_g$, the irregularity of the boundary starts to affect the statistics of
<figure><img src="image_11.png" /><figcaption>Figure 10: Comparison of empirical PDFs (in log-log scale) of the boundary local time $\ell_t$ normalized by its mean value, obtained in the Koch snowflake with $L = 2$, $D = 1$, and $N = 10^8$ particles by Monte Carlo simulations at different times and generations: $g = 10$, $t = 10^{-7}$ (blue line), $g = 9$, $t = 9 \times 10^{-7}$ (orange dot), $g = 8$, $t = 9^2 \times 10^{-7}$, and $g = 7$, $t = 9^3 \times 10^{-7}$ (red diamond).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Mean boundary local time (rescaled by $\sqrt{Dt}$) obtained in the Koch snowflake with $L = 2$, $D = 1$, and $N = 10^8$ particles by Monte Carlo simulations for different times $t$ and different generations $g$. The dashed line presents an empirical fitting.</figcaption></figure>

encounters, yielding an increase of the rescaled mean value. This increase is observed up to $t/t_g \approx 10^4$. Fitting $\langle \ell_t \rangle_0 / \sqrt{Dt} = f(t/t_g)$ in this region, we find an approximation

$$f(z) = \frac{6}{\sqrt{\pi}} \left[ 1 + \frac{\ln(0.15z)}{6} \Theta(z - 1/0.15) \right], \quad (31)$$
<figure><img src="image_13.png" /><figcaption>Figure 12: Comparison of local persistence exponents $\alpha_q(t|x_0)$ obtained in the sixth generation $\Omega_6$ of the Koch snowflake with $L = 2$, $D = 1$, and $N = 10^8$ particles by Monte Carlo simulations for different reactivities: $q = \infty$ (thin black line), $q = 10^5$ (blue dots), $q = 10^4$ (orange plus), $q = 10^3$ (green crosses), $q = 10^2$ (red up triangles), $q = 10^1$ (violet down triangles), and $q = 10^0$ (brown squares). The median time $T_q$ for each reactivity (from $q = 10^5$ to $q = 10^0$) is indicated by the vertical line with the corresponding color: $1.09 \times 10^{-7}$, $1.44 \times 10^{-7}$, $5.06 \times 10^{-7}$, $8.66 \times 10^{-6}$, $3.13 \times 10^{-4}$, $8.60 \times 10^{-3}$.</figcaption></figure>

where $\Theta(x)$ is the Heaviside step function. When $t/t_g$ exceeds $10^4$, the rescaled mean value starts to decrease. This behavior qualitatively agrees with what we observed for the simplified case of a wedge with channels: the mean boundary local time first increases and then decreases (Fig. 8). Remarkably, the effect of the Koch snowflake complexity is rather weak: the ratio $\langle l_t \rangle_0 / \sqrt{Dt}$ varies by a factor of 2, whereas $t/t_g$ changes over many orders of magnitude. A similar behavior was observed for other Koch snowflakes of different angles $\alpha$ (results are not shown). As $\alpha$ increases, the maximum is shifted to larger times, e.g., $t/t_g \approx 10^6$ for $\alpha = 3\pi/4$.

### 4.3. Local persistence exponents

In many complex systems, the survival probability exhibits a power-law decay (see [81–83] and references therein). In our previous work [59], we studied the long-time asymptotic behavior of the survival probability $S_\infty(t|x_0)$ in the Koch snowflake with a perfectly reactive boundary. For this purpose, the local persistence exponent (LPE) was introduced as the negative logarithmic derivative of the survival probability that we extend here to the case of finite reactivity:

$$ \alpha_q(t|x_0) = - \frac{\partial \ln S_q(t|x_0)}{\partial \ln t}. \qquad (32) $$
Due to the self-similarity of the Koch snowflake boundary, $\alpha_{\infty}(t|x_0)$ was shown to exhibit log-periodic oscillations in time. Here we inspect how the finite reactivity of the boundary can affect this behavior.

We compute the LPE of the survival probability $S_q(t|x_0)$ by estimating the latter via the EFS approach. The starting point is located at $x_0 = (x_0, y_0 + 10^{-3})$, where $(x_0, y_0)$ is the bottom vertex of the Koch snowflake as shown in Fig. 1b. Here, the starting point is not located at the vertex, in order to enable comparison with $\alpha_{\infty}(t|x_0)$ for the perfectly reactive case. To avoid too long computational time for small values of $q$, we focus on the generation $g = 6$. Figure 12 compares LPEs for different reactivities $q$ from infinity to 1.

As discussed in [59], $\alpha_{\infty}(t|x_0)$ exhibits a transient regime of monotonous growth at very short times, log-periodic oscillations at intermediate times due to self-similarity of the boundary, and further linear increase with $t$ at long times due to confinement in a bounded domain. We observe here similar trends for $\alpha_q(t|x_0)$. When the reactivity is very high ($q = 10^5$), the particle reacts after first few encounters with the boundary, so that $S_q(t|x_0)$ is close to $S_{\infty}(t|x_0)$, implying $\alpha_q(t|x_0) \approx \alpha_{\infty}(t|x_0)$. A decrease of reactivity extends the transient regime to larger and larger times. Moreover, the ultimate linear increase, $\alpha_q(t|x_0) \propto t$ starts at $t \approx L^2/D$, independently of $q$. As a consequence, the range of log-periodic oscillations is reduced as $q$ decreases. For instance, there is no such an intermediate regime for $q = 1$. Note that a natural timescale of the transient regime is the median reaction time $T_q$, at which half of the particles survive $S_q(T_q|x_0) = 1/2$. In fact, a significant decay of the survival probability is expected at $t \gg T_q$.

## 5. Discussion and conclusion

In this paper, we studied the boundary local time in polygonal domains. First, our theoretical framework outlined the statistics of the boundary local time in a wedge, including the mean value, the variance, and the asymptotic behavior of the probability density. In fact, Eqs. (10, 11) give the mean value and the variance of the boundary local time when the starting point is located at the vertex, whereas their extensions to arbitrary starting points are given in Appendix C; note that the prefactors $v_{\alpha}$ in Eq. (11) were found explicitly for the angles $\alpha = \pi/n$ with an integer $n$, whereas their computation for any $\alpha$ is possible via coupled PDEs between moments $\langle \ell_t^k \rangle_{x_0}$ from Appendix A. The asymptotic behavior (22) of the probability density $\rho_{\alpha}(\ell, t|x_0)$ in the limit $\ell \to \infty$ was derived for wedges of angle $\alpha \le \pi$ from the duality between the Steklov problem and the Robin Laplacian problem. As the small-$\ell$ behavior relies on all Steklov eigenmodes, we could not establish this behavior rigorously, but we provided its conjectural form.

To deal with polygonal domains, we developed an efficient escape-from-a-sector approach, which relies on the approximate solution of the escape problem for a sector. Comparison of simulated values of $\ell_t$ with theoretical predictions in wedges helped us to validate the EFS approach and to access its accuracy. The EFS approach was also
verified by computing the spread harmonic measure distribution in the Koch snowflake and comparing with an alternative computation by a FEM.

The non-monotonous behavior of the mean boundary local time in complex domains presents one of the main numerical results of the paper. A compromise between an increase of the perimeter and its reduced accessibility due to larger spaces for diffusion was identified as its origin. Indeed, the larger perimeter tends to increase the mean boundary local time, while the larger space for diffusion diminishes the chance of encounters on the surface in a given time. Such a non-monotonous dependence was first observed for a simple model of the wedge of angle $\pi/2$ with two rectangular channels and then for finite generations of the Koch snowflake. In the latter case, a clear increase of the mean boundary local time was noticed for $t/t_g \in (0, 10^4)$, which is followed by a decrease. The statistics of encounters is much richer due to the self-similarity of the boundary. When the starting point is located at the bottom vertex, changing the generation $g$ is equivalent to rescaling time $t$, yielding identical PDFs of $\ell_t$ after normalization by the mean value $\langle \ell_t \rangle_0$.

Self-similarity of the Koch snowflake boundary leads to log-periodic oscillations of the local persistence exponent of the survival probability. We examined how the finite reactivity affects this behavior. Three regimes were distinguished: (i) a transient growth at very short times; (ii) log-oscillations at intermediate times; (iii) linear growth at long times due to confinement in a bounded domain. We introduced the median time $T_q$ as a suitable timescale for the transient regime. When $t \lesssim T_q$, the number of encounters of the diffusing particle with the boundary is not enough to ensure the reaction event, and the survival probability is close to 1, yielding $\alpha_q(t|x_0) \approx 0$. In turn, the asymptotic decay of the survival probability is observed at times $t \gg T_q$. As a consequence, when the reactivity $q$ is small, $T_q$ is large, and the intermediate regime with log-periodic oscillations disappears. In practice, both geometry and reactivity affect the survival probability, its asymptotic behavior, and the median time $T_q$. Their systematic analysis presents an interesting perspective of this work.

## Acknowledgements

D.S.G. acknowledges the Simons Foundation for supporting his sabbatical sojourn in 2024 at the CRM (CNRS - University of Montreal, Canada), as well as the Alexander von Humboldt Foundation for support within a Bessel Prize award. Y.Y acknowledges Adrien Chaigneau and Xiaozhen Wang for fruitful discussions about the asymptotic behavior. We also thank Adrien Chaigneau for his computations of the spread harmonic measure on the Koch snowflake by a FEM (Fig. 6) for validation purposes.

## Appendix A. Coupled PDEs for the moments

The moments of the boundary local time can be found from the propagator $G_0(x, t|x_0)$ via a hierarchical set of PDEs. Their derivation relies on the relation (4) between the
survival probability and the probability density function of $\ell_t$. The survival probability $S_q(t|x_0)$ that satisfies Eqs. (5), admits the following Taylor expansion

$$S_q(t|x_0) = \sum_{k=0}^{\infty} \frac{(-q)^k L_k(t|x_0)}{k!}, \qquad (A.1)$$

where $L_k(t|x_0) = \langle \ell_t^k \rangle_{x_0}$ denotes the k-th order moment of $\ell_t$. Substituting this Taylor expansion into Eqs. (5), we get a sequence of PDEs for the moments with $k = 1, 2, \dots$

$$\partial_t L_k(t|x_0) = D \Delta L_k(t|x_0) \quad \text{in } \Omega, \quad \partial_n L_k(t|x_0) = k L_{k-1}(t|x_0) \quad \text{on } \partial\Omega, \quad L_k(0|x_0) = 0. \qquad (A.2)$$

Note that for $k = 0, L_0(t|x_0) = 1$. Equivalently, the Laplace transform of $L_k(t|x_0)$, $\tilde{L}_k(p|x_0) = \int_0^\infty dt e^{-pt} L_k(t|x_0)$, satisfies

$$(p - D\Delta)\tilde{L}_k(p|x_0) = 0 \quad \text{in } \Omega, \quad \partial_n \tilde{L}_k(p|x_0) = k \tilde{L}_{k-1}(p|x_0) \quad \text{on } \partial\Omega. \qquad (A.3)$$

Let $\tilde{G}_0(x, p|x_0)$ denote the Laplace transform of the propagator $G_0(x, t|x_0)$ with the Neumann boundary condition that satisfies

$$(p - D\Delta)\tilde{G}_0(x, p|x_0) = \delta(x - x_0) \quad \text{in } \Omega, \quad \partial_n \tilde{G}_0(x, p|x_0) = 0 \quad \text{on } \partial\Omega. \qquad (A.4)$$

Multiplying Eq. (A.3) by $\tilde{G}_0(x, p|x_0)$, multiplying Eq. (A.4) by $\tilde{L}_k(p|x_0)$, subtracting these equations, integrating them over $x_0 \in \Omega$, and using the Green's formula, we get

$$\tilde{L}_k(p|x) = k \int_{\partial\Omega} d\mathbf{x}_0 D\tilde{G}_0(\mathbf{x}, p|\mathbf{x}_0) \tilde{L}_{k-1}(p|\mathbf{x}_0). \qquad (A.5)$$

In the time domain, this equation implies:

$$L_k(t|x) = k \int_{\partial\Omega} d\mathbf{x}_0 \int_0^t dt' D G_0(\mathbf{x}, t'|\mathbf{x}_0) L_{k-1}(t-t'|\mathbf{x}_0). \qquad (A.6)$$

As $L_0(t|x_0) = 1$ and thus $\tilde{L}_0(p|x_0) = 1/p$, we have

$$\tilde{L}_1(p|x) = \frac{1}{p} \int_{\partial\Omega} d\mathbf{x} \, D\tilde{G}_0(\mathbf{x}, p|x_0), \qquad (A.7)$$

and we therefore retrieve the relation (14):

$$L_1(t|x_0) = \int_0^t dt' \int_{\partial\Omega} d\mathbf{x} \, D G_0(\mathbf{x}, t'|x_0). \qquad (A.8)$$

In turn, the second moment reads

$$L_2(t|x_0) = 2 \int_0^t dt' \int_{\partial\Omega} d\mathbf{x} \, D G_0(\mathbf{x}, t'|x_0) L_1(t-t'|\mathbf{x}) \qquad (A.9)$$

$$= 2D^2 \int_{\partial\Omega} d\mathbf{x}_1 \int_{\partial\Omega} d\mathbf{x}_2 \int_0^t dt_1 \int_0^{t_1} dt_2 G_0(\mathbf{x}_2, t_2|x_1) G_0(\mathbf{x}_1, t-t_1|x_0). \qquad (A.10)$$
If the boundary $\partial\Omega$ of the domain is bounded, one can use the spectral decomposition of $\tilde{G}_0(\boldsymbol{x}, p|\boldsymbol{x}_0)$ over the Steklov eigenfunctions $V_k^{(p)}$, which reads for $\boldsymbol{x} \in \partial\Omega$ as [7]

$$D\tilde{G}_0(\boldsymbol{x}, p|\boldsymbol{x}_0) = \sum_k \frac{V_k^{(p)}(\boldsymbol{x}_0)V_k^{(p)}(\boldsymbol{x})}{\mu_k^{(p)}}. \qquad (A.11)$$

The orthogonality of Steklov eigenfunctions on the boundary allows one to compute integrals over $\partial\Omega$ in Eq. (A.5) and thus to express the Laplace transform of the k-th moment as

$$\tilde{L}_k(p|\boldsymbol{x}_0) = \frac{1}{p} \sum_k \frac{V_k^{(p)}(\boldsymbol{x}_0)}{[\mu_k^{(p)}]^k} \int_{\partial\Omega} d\boldsymbol{x} V_k^{(p)}(\boldsymbol{x}), \qquad (A.12)$$

i.e., we retrieved the spectral expansion from [60]. However, this computation is not applicable when the boundary $\partial\Omega$ is unbounded (as in the case of wedges), and one needs to evaluate the integrals in Eq. A.6.

## Appendix B. Explicit results for a quadrant

For completeness, we provide some explicit results for the quadrant $\mathbb{R}_+^2$, i.e., the wedge of angle $\pi/2$. The following computation is elementary and relies on the independence of horizontal and vertical displacements that reduces the analysis to one-dimensional problems. Despite their simplicity, we are not aware of earlier references presenting these results.

We recall that the probability density of the boundary local time $\ell_t$ on the positive semi-axis reads as

$$\rho_{\pi}(\ell, t|x_0) = \text{erf}\left(\frac{x_0}{\sqrt{4Dt}}\right) \delta(\ell) + \frac{e^{-(\ell+x_0)^2/(4Dt)}}{\sqrt{\pi Dt}}, \qquad (B.1)$$

where the first singular term, $\delta(\ell)$, accounts for trajectories that do not hit the boundary up to time $t$, with $\text{erf}(x_0/\sqrt{4Dt})$ being the associated survival probability. As a consequence, we have

$$\begin{aligned} \rho_{\pi/2}(\ell, t|x_0, y_0) &= \int_0^\infty d\ell_1 \rho_\pi(\ell_1, t|x_0) \int_0^\infty d\ell_2 \rho_\pi(\ell_2, t|y_0) \delta(\ell_1 + \ell_2 - \ell) \\ &= S_x S_y \delta(\ell) + \frac{S_x}{\sqrt{\pi Dt}} e^{-(\ell+y_0)^2/(4Dt)} + \frac{S_y}{\sqrt{\pi Dt}} e^{-(\ell+x_0)^2/(4Dt)} \\ &\quad + \frac{2}{\pi \sqrt{Dt}} \int_0^\infty dz_1 \int_0^\infty dz_2 e^{-(z_1+\xi_1)^2-(z_2+\xi_2)^2} \delta(z_1 + z_2 - \ell/\sqrt{4Dt}), \end{aligned}$$

where $\xi_1 = x_0/\sqrt{4Dt}$, $\xi_2 = y_0/\sqrt{4Dt}$, $S_x = \text{erf}(x_0/\sqrt{4Dt})$, and $S_y = \text{erf}(y_0/\sqrt{4Dt})$.
Evaluating the last integral, we have

$$ \begin{aligned} \rho_{\pi/2}(\ell, t|x_0, y_0) = & \operatorname{erf}\left(\frac{x_0}{\sqrt{4Dt}}\right) \operatorname{erf}\left(\frac{y_0}{\sqrt{4Dt}}\right) \delta(\ell) \\ & + \frac{e^{-(\ell+y_0)^2/(4Dt)}}{\sqrt{\pi Dt}} \operatorname{erf}\left(\frac{x_0}{\sqrt{4Dt}}\right) + \frac{e^{-(\ell+x_0)^2/(4Dt)}}{\sqrt{\pi Dt}} \operatorname{erf}\left(\frac{y_0}{\sqrt{4Dt}}\right) \\ & + \frac{e^{-(\ell+x_0+y_0)^2/(8Dt)}}{\sqrt{2\pi Dt}} \left[ \operatorname{erf}\left(\frac{\ell-y_0+x_0}{\sqrt{8Dt}}\right) + \operatorname{erf}\left(\frac{\ell+y_0-x_0}{\sqrt{8Dt}}\right) \right]. \end{aligned} \quad (B.2) $$

As the boundary local times on the horizontal and vertical axes are independent, we get

$$ \langle \ell_t \rangle_{x_0} = \frac{2\sqrt{Dt}}{\sqrt{\pi}} \left[ e^{-x_0^2/(4Dt)} - \frac{\sqrt{\pi} x_0}{\sqrt{4Dt}} \operatorname{erfc}\left(\frac{x_0}{\sqrt{4Dt}}\right) + e^{-y_0^2/(4Dt)} - \frac{\sqrt{\pi} y_0}{\sqrt{4Dt}} \operatorname{erfc}\left(\frac{y_0}{\sqrt{4Dt}}\right) \right]. \quad (B.3) $$

This relation agrees with our general expression (C.21), with $m = 2$, in which two terms in the sum can be associated to $y_0 = r_0 \sin \theta_0$ and $x_0 = r_0 \sin \theta_1 = r_0 \cos \theta_0$.

We also note that the probability density $\rho_\alpha(\ell, t|x_0)$ of the boundary local time $\ell_t$ is tightly related to the probability density $U_\alpha(\ell, t|x_0)$ of the first-crossing time $\mathcal{T}_\ell = \inf\{t > 0 : \ell_t > \ell\}$ of a threshold $\ell$ [7]:

$$ \int_{\ell}^{\infty} d\ell' \rho_{\alpha}(\ell', t|x_0) = \mathbb{P}_{x_0}\{\ell_t > \ell\} = \mathbb{P}_{x_0}\{\mathcal{T}_\ell < t\} = \int_{0}^{t} dt' U_\alpha(\ell, t'|x_0), \quad (B.4) $$

i.e.,

$$ U_{\alpha}(\ell, t|x_0) = \int_{\ell}^{\infty} d\ell' \partial_{t} \rho_{\alpha}(\ell', t|x_0). \quad (B.5) $$

For instance, for the wedge of angle $\pi$, one has

$$ U_{\pi}(\ell, t|0) = \frac{\ell e^{-\ell^2/(4Dt)}}{\sqrt{4\pi Dt^3}}. \quad (B.6) $$

Similarly, we get for the quadrant

$$ U_{\pi/2}(\ell, t|0) = \frac{\ell e^{-\ell^2/(8Dt)} \operatorname{erf}(\ell/\sqrt{8Dt})}{\sqrt{2\pi Dt^3}}. \quad (B.7) $$

In both cases, the mean first-crossing time is infinite, as expected from the divergence of the mean first-passage time.

## Appendix C. Mean boundary local time and its variance in wedges

In this Appendix, we elaborate on the computation of the mean boundary local time (Appendix C.2) and its variance (Appendix C.3). These computations require knowledge of the propagator for the wedge. Even though the computation of the propagator is standard, we reproduce it in Appendix C.1 for completeness.
Appendix C.1. Propagator in the wedge

Let us compute the propagator $G_0(x, t|x_0)$ in a wedge of angle $\alpha$ with Neumann boundary condition. Its Laplace transform, $\tilde{G}_0(x, p|x_0) = \int_0^\infty dt e^{-pt} G_0(x, t|x_0)$, satisfies

$$ (p - D\Delta)\tilde{G}_0(x, p|x_0) = \delta(x - x_0) = r_0^{-1}\delta(r - r_0)\delta(\theta - \theta_0) \quad \text{in } \Omega, \\ \partial_n\tilde{G}_0(x, p|x_0) = 0 \quad \text{on } \partial\Omega. \qquad (C.1) $$

One can search for its solution in the form that respects the boundary condition:

$$ \tilde{G}_0(x, p|x_0) = \frac{1}{\alpha}\tilde{g}_0(r, p|r_0) + \frac{2}{\alpha}\sum_{n=1}^{\infty}\cos(\nu_n\theta)\cos(\nu_n\theta_0)\tilde{g}_n(r, p|r_0), \qquad (C.2) $$

with $\nu_n = \pi n / \alpha$ and unknown radial functions $\tilde{g}_n(r, p|r_0)$. Substitution of this form into the above equation yields

$$ \mathcal{L}_0 \tilde{g}_0(r, p|r_0) + 2 \sum_{n=1}^{\infty} \cos(\nu_n \theta) \cos(\nu_n \theta_0) \mathcal{L}_n \tilde{g}_n(r, p|r_0) = \frac{\alpha}{Dr_0} \delta(r - r_0) \delta(\theta - \theta_0), \qquad (C.3) $$

where $\mathcal{L}_n = p/D - (\partial_r^2 + r^{-1}\partial_r - \nu_n^2 r^{-2})$. This equation can be satisfied by imposing

$$ \mathcal{L}_n \tilde{g}_n(r, p|r_0) = \frac{1}{Dr_0} \delta(r - r_0). \qquad (C.4) $$

For each $n$, the solution of this equation can be found by solving the homogeneous equation for $r < r_0$ and $r > r_0$ and then matching two solutions. One gets

$$ \tilde{g}_n(r, p|r_0) = \frac{1}{D} I_{\nu_n}(r_< \sqrt{p/D}) K_{\nu_n}(r_> \sqrt{p/D}), \qquad (C.5) $$

where $r_< = \min\{r, r_0\}$, $r_> = \max\{r, r_0\}$, and $I_{\nu_n}(z)$ and $K_{\nu_n}(z)$ are the modified Bessel functions of the first and second kind. The inverse Laplace transform of this function yields [84]

$$ g_n(r, t|r_0) = \frac{1}{2Dt} e^{-(r^2+r_0^2)/(4Dt)} I_{\nu_n}(rr_0/(2Dt)). \qquad (C.6) $$

We conclude that

$$ G_0(x, t|x_0) = \frac{g_0(r, t|r_0)}{\alpha} + \frac{2}{\alpha} \sum_{n=1}^{\infty} \cos(\nu_n \theta) \cos(\nu_n \theta_0) g_n(r, t|r_0). \qquad (C.7) $$

Appendix C.2. Mean boundary local time

The substitution of Eq. (C.7) into Eq. (14) yields the mean boundary local time

$$ \begin{aligned} \langle \ell_t \rangle_{\mathbf{x}_0} &= D \int_0^t dt' \int_0^\infty dr \left[ G_0(\mathbf{x}, t'|\mathbf{x}_0)|_{\theta=0} + G_0(\mathbf{x}, t'|\mathbf{x}_0)|_{\theta=\alpha} \right] \\ &= \frac{D}{\alpha} \sum_{n=0}^{\infty} \epsilon_n \cos(\nu_n \theta_0) (1 + (-1)^n) \int_0^t dt' e^{-r_0^2/(8Dt')} \frac{\sqrt{\pi}}{2\sqrt{4Dt'}} I_{\nu_n/2}(r_0^2/(8Dt')), \end{aligned} \qquad (C.8) $$
with $\epsilon_n = 2 - \delta_{n,0}$, and we used

$$ \int_0^\infty dz \, e^{-z^2} I_{\nu_n}(az) = \frac{\sqrt{\pi}}{2} e^{a^2/8} I_{\nu_n/2}(a^2/8). \qquad (C.9) $$

At $r_0 = 0$, only the term with $n = 0$ contributes, yielding

$$ \frac{1}{\alpha} \int_{0}^{\infty} dr \, 2g_{0}(r, t|r_{0}) = \frac{\sqrt{\pi}}{\alpha \sqrt{Dt}}, \qquad (C.10) $$

from which the integral over $t$ implies

$$ \langle l_t \rangle_0 = \frac{\pi}{\alpha} \frac{2\sqrt{Dt}}{\sqrt{\pi}}. \qquad (C.11) $$

We therefore retrieved the expression obtained in the main text by probabilistic and symmetry arguments.

In turn, for $r_0 > 0$, one can rewrite the above expression as

$$ \langle l_t \rangle_{x_0} = \frac{r_0 \sqrt{\pi}}{2\sqrt{2\alpha}} \sum_{n=0}^{\infty} \epsilon_n \cos(2\nu_n \theta_0) \int_{r_0^2/(8Dt)}^{\infty} \frac{dz}{z^{3/2}} e^{-z} I_{\nu_n}(z). \qquad (C.12) $$

Note that an accurate numerical computation of this sum may require taking a significant number of terms.

*Special case.* In the special case $\alpha = \pi/m$ with an integer $m = 1, 2, \dots$, one can further simplify Eq. (C.12) to get a closed-form expression. For this purpose, we use the following representation of the modified Bessel function of the first kind:

$$ I_{\nu_n}(z) = \frac{1}{\pi} \int_0^\pi d\theta \cos(\nu_n \theta) e^{z \cos \theta} - \frac{\sin(\nu_n \pi)}{\pi} \int_0^\infty dx \, e^{-z \cosh x - \nu_n x}. \qquad (C.13) $$

Since $\nu_n = \pi n/\alpha = mn$ is integer, the second term vanishes, whereas the substitution of the first term into Eq. (C.12) yields

$$ \langle l_t \rangle_{x_0} = \frac{2\sqrt{Dt}}{\sqrt{\pi}} L^{(1)}(\eta), \qquad (C.14) $$

where $\eta = r_0^2/(8Dt)$, and

$$ L^{(1)}(\eta) = \frac{\sqrt{\eta}}{2\alpha} \int_0^\pi d\theta \int_\eta^\infty \frac{dz}{z^{3/2}} e^{-z} e^{z \cos \theta} \sum_{n=0}^\infty \epsilon_n \cos(2\nu_n \theta_0) \cos(\nu_n \theta). \qquad (C.15) $$
Using the identity

$$ \int_{x}^{\infty} \frac{dz}{z^{3/2}} e^{-z} = \frac{2e^{-x}}{\sqrt{x}} - \sqrt{\pi} \operatorname{erfc}(\sqrt{x}), \qquad (C.16) $$

we get

$$ L^{(1)}(\eta) = \frac{1}{\alpha} \int_{0}^{\pi} d\theta \left[ e^{-\eta(1-\cos\theta)} - \sqrt{\pi} \sqrt{\eta(1-\cos\theta)} \operatorname{erfc}(\sqrt{\eta(1-\cos\theta)}) \right] \times \sum_{n=0}^{\infty} \epsilon_n \cos(2\pi n \theta_0 / \alpha) \cos(\pi n \theta / \alpha). \quad (C.17) $$

Due to the reflection symmetry with respect to the central ray of the wedge, we can assume that $2\theta_0 \le \alpha$. Setting $\xi_0 = 2\theta_0/\alpha$, we see that $0 \le \xi_0 \le 1$. In turn, $\xi = \theta/\alpha$ can take any value from 0 to $m$ given that $\theta$ varies up to $\pi$ in Eq. (C.17). We aim at applying the following identity

$$ \sum_{n=0}^{\infty} \epsilon_n \cos(\pi n \xi_0) \cos(\pi n \xi) = \delta(\xi_0 - \xi) \quad (0 < \xi_0, \xi < 1). \qquad (C.18) $$

For this purpose, the integral over $\theta$ from 0 to $\pi$ can be split into $m$ integrals over the intervals $(\alpha j, \alpha(j+1))$, with $j = 0, 1, \dots, m-1$. To apply the above identity, we introduce the angles

$$ \theta_j = \begin{cases} \theta_0 + \frac{1}{2}\alpha j & \text{if } j \text{ is even,} \\ \frac{1}{2}\alpha(j+1) - \theta_0 & \text{if } j \text{ is odd,} \end{cases} \qquad (C.19) $$

for $j = 0, 1, 2, \dots, m-1$. These angles determine the mirror reflections of the starting point $(r_0, \theta_0)$ with respect to the reflecting boundary of the wedge. Applying Eq. (C.18) with a suitable shift of $\xi$ to ensure that $\xi - \xi_j \in (0, 1)$ for each interval $(\alpha j, \alpha(j+1))$, we get

$$ L^{(1)}(\eta) = \sum_{j=0}^{m-1} \left[ e^{-\eta(1-\cos(2\theta_j))} - \sqrt{\pi} \sqrt{\eta(1-\cos(2\theta_j))} \operatorname{erfc}\left(\sqrt{\eta(1-\cos(2\theta_j))}\right) \right], \quad (C.20) $$

from which

$$ \langle l_t \rangle_{x_0} = \frac{2\sqrt{Dt}}{\sqrt{\pi}} \sum_{j=0}^{m-1} \left[ e^{-r_0^2 \sin^2 \theta_j / (4Dt)} - \sqrt{\pi} \frac{r_0 \sin \theta_j}{\sqrt{4Dt}} \operatorname{erfc} \left( \frac{r_0 \sin \theta_j}{\sqrt{4Dt}} \right) \right]. \quad (C.21) $$

In the limit $r_0 \to 0$, the sum is simply equal to $m = \pi/\alpha$, and we retrieve Eq. (C.11).

**General case.** Let us first focus on the case $\alpha > \pi$. As previously, we assume that $\theta_0 \le \alpha/2$ (otherwise $\theta_0$ can be replaced by $\alpha - \theta_0$ by symmetry). As $\pi/\alpha$ is not an integer anymore, both integrals in Eq. (C.13) have to be considered. The contribution of the first integral is still given by Eq. (C.20) with $m = 1$ so that only one term is
present. However, the contribution of this single term is canceled if $2\theta_0 > \pi$, due to the Dirac distribution in Eq. (C.18) and the integration over $\theta$ from 0 to $\pi$. We conclude that

$$L^{(1)}(\eta) = \Theta(\pi - 2\theta_0) \left[ e^{-2\eta \sin^2 \theta_0} - \sqrt{2\pi\eta} \sin \theta_0 \operatorname{erfc}(\sqrt{2\eta} \sin \theta_0) \right]. \quad (C.22)$$

For the second term, we have

$$L^{(2)}(\eta) = \frac{\sqrt{\eta}}{2\alpha} \int_0^\infty dx \int_\eta^\infty \frac{dz}{z^{3/2}} e^{-z} e^{-z \cosh x} \sum_{n=1}^\infty 2 \cos(2\nu_n \theta_0) \sin(\nu_n \pi) e^{-\nu_n x}.$$

Using the identity

$$\sum_{n=1}^{\infty} \sin(n\beta_{\pm})e^{-n\delta} = \frac{e^{-\delta} \sin \beta_{\pm}}{1 - 2 \cos \beta_{\pm} e^{-\delta} + e^{-2\delta}}, \quad (C.23)$$

with $\delta = \pi x/\alpha$ and $\beta_{\pm} = \pi(\pi \pm 2\theta_0)/\alpha$, we find

$$L^{(2)}(\eta) = \frac{1}{\alpha} \int_{0}^{\infty} dx \left[ \frac{e^{-\pi x/\alpha} \sin \beta_{+}}{1 - 2 \cos \beta_{+} e^{-\pi x/\alpha} + e^{-2\pi x/\alpha}} + \frac{e^{-\pi x/\alpha} \sin \beta_{-}}{1 - 2 \cos \beta_{-} e^{-\pi x/\alpha} + e^{-2\pi x/\alpha}} \right] \\ \times \left[ e^{-\eta(1+\cosh x)} - \sqrt{\pi\eta(1+\cosh x)} \operatorname{erfc}(\sqrt{\eta(1+\cosh x)}) \right].$$

Changing the integration variable, $z = e^{-\pi x/\alpha}$, we have

$$L^{(2)}(\eta) = \frac{1}{\pi} \int_{0}^{1} dz \left( \frac{\sin \beta_{+}}{1 - 2z \cos \beta_{+} + z^2} + \frac{\sin \beta_{-}}{1 - 2z \cos \beta_{-} + z^2} \right) \\ \times \left[ e^{-(z^{\gamma} + z^{-\gamma})^2 \eta/2} - \sqrt{\pi\eta/2}(z^{\gamma} + z^{-\gamma}) \operatorname{erfc}((z^{\gamma} + z^{-\gamma})\sqrt{\eta/2}) \right],$$

where $\gamma = \alpha/(2\pi)$. Combining this expression with Eq. (C.22), we get a relatively simple expression for an accurate computation of the mean boundary local time:

$$\langle l_t \rangle_{x_0} = \frac{2\sqrt{Dt}}{\sqrt{\pi}} (L^{(1)}(\eta) - L^{(2)}(\eta)). \quad (C.24)$$

Finally, when $\alpha < \pi$ but the ratio $\pi/\alpha$ is not integer, the second contribution $L^{(2)}(\eta)$ remains unchanged. In turn, the first contribution admits a form similar to Eq. (C.20), in which $m = 1 + \lfloor \pi/\alpha^{+0} \rfloor$, where $\lfloor x \rfloor$ is the integer part of $x$ (i.e., the greatest integer less than or equal to $x$); here $\alpha^{+0}$ is the convention to take the limit of $\lfloor \pi/(\alpha + \epsilon) \rfloor$ as $\epsilon \to 0$ to naturally incorporate the case when $\alpha = \pi/m$ with an integer $m$. The only difference with Eq. (C.20) is that if $2\theta_j$ exceeds $\pi$, the contribution of the corresponding term should be canceled, as in Eq. (C.22). We conclude that

$$L^{(1)}(\eta) = \sum_{j=0}^{m-1} \Theta(\pi/2 - \theta_j) \left[ e^{-2\eta \sin^2 \theta_j} - \sqrt{2\pi\eta} \sin \theta_j \operatorname{erfc}(\sqrt{2\eta} \sin \theta_j) \right], \quad (C.25)$$

and Eq. (C.24) is now valid for any $0 < \alpha \le 2\pi$. We recall that $\theta_0$ should be understood as $\min\{\theta_0, \alpha - \theta_0\}$.
Appendix C.3. Computation of the variance for $\alpha = \pi/m$

Since the first moment $L_1(t|x_0) = \langle l_t \rangle_{x_0}$ is known via Eq. (C.12), we can also compute the second moment by evaluating the integrals in Eq. (A.9). However, this computation in the general case is tedious. For this reason, we focus on the starting point $x_0 = 0$ and consider only the special case $\alpha = \pi/m$ with an integer $m$, for which the first moment admits a simple form (C.21). In this case, Eq. (A.9) reads

$$L_2(t|0) = 2D \int_0^t dt' \int_0^\infty dr \left[ G_0(r, 0, t'|0) L_1(t-t'|r, 0) + G_0(r, \alpha, t'|0) L_1(t-t'|r, \alpha) \right], \quad (C.26)$$

where two terms represent the integrals over two rays at $\theta = 0$ and $\theta = \alpha$, written in polar coordinates. At the starting point $x_0 = 0$, the propagator has a simple form, given by the first term in Eq. (C.7). As both the propagator and the first moment are symmetric with respect to the central ray of the wedge, the two terms in Eq. (C.12) are identical. Substituting Eq. (C.6) and Eq. (C.21), we get thus

$$L_2(t|0) = 4D \int_0^t dt' \int_0^\infty dr \frac{e^{-r^2/(4Dt')}}{2Dt'\alpha} \left\{ \frac{2\sqrt{D(t-t')}}{\sqrt{\pi}} \sum_{j=0}^{m-1} \left[ e^{-r^2 \sin^2 \theta_j / (4D(t-t'))} \right. \right. \\ \left. \left. - \frac{\sqrt{\pi}r \sin \theta_j}{\sqrt{4D(t-t')}} \operatorname{erfc} \left( \frac{r \sin \theta_j}{\sqrt{4D(t-t')}} \right) \right] \right\}, \quad (C.27)$$

where $\theta_j$ were defined in Eq. (C.19) with $\theta_0 = 0$. Changing the integration variables $z = r/\sqrt{4Dt'}$ and $\tau = t'/t$, we have

$$L_2(t|0) = \frac{8Dt}{\sqrt{\pi}\alpha} \int_0^1 \frac{d\tau}{\sqrt{\tau}} \sqrt{1-\tau} \int_0^\infty dz e^{-z^2} \sum_{j=0}^{m-1} \left[ e^{-z^2 \sin^2 \theta_j \tau / (1-\tau)} \right. \\ \left. - z \sin \theta_j \sqrt{\frac{\pi\tau}{1-\tau}} \operatorname{erfc} \left( z \sin \theta_j \sqrt{\frac{\tau}{1-\tau}} \right) \right]. \quad (C.28)$$

Evaluation of the integrals over $z$ yields

$$L_2(t|0) = \frac{4Dt}{\alpha} \sum_{j=0}^{m-1} \int_0^1 \frac{d\tau}{\sqrt{\tau}} \sqrt{1-\tau} \left\{ \sqrt{\frac{1-\tau}{1-\tau \cos^2 \theta_j}} - \sin \theta_j \sqrt{\frac{\tau}{1-\tau}} \left( 1 - \frac{A_\tau}{\sqrt{1+A_\tau^2}} \right) \right\}, \quad (C.29)$$

where $A_\tau = \sin \theta_j \sqrt{\tau/(1-\tau)}$. We get then

$$L_2(t|0) = \frac{4Dt}{\alpha} \sum_{j=0}^{m-1} \int_0^1 \frac{d\tau}{\sqrt{\tau}} \left\{ \frac{1-\tau}{\sqrt{1-\tau \cos^2 \theta_j}} - \sin \theta_j \sqrt{\tau} \left( 1 - \frac{\sqrt{\tau} \sin \theta_j}{\sqrt{1-\tau \cos^2 \theta_j}} \right) \right\} \\ = \frac{4Dt}{\alpha} \sum_{j=0}^{m-1} \left\{ -\sin \theta_j + \int_0^1 \frac{d\tau}{\sqrt{\tau}} \sqrt{1-\tau \cos^2 \theta_j} \right\} = \frac{4Dt}{\alpha} \sum_{j=0}^{m-1} \frac{\pi/2 - \theta_j}{\cos \theta_j} \quad (C.30)$$
(note that if $\theta_j = \pi/2$, the ratio is set to 1). Recalling Eq. (C.19), we obtain the explicit form of the second moment

$$L_2(t|0) = 4Dt \cdot \sigma_m, \qquad (C.31)$$

with the prefactor $\sigma_m$ that depends on the parity of $m$:

$$\sigma_m = m \left[ \frac{1}{2} + \frac{1}{\pi} + \sum_{j=1}^{(m-2)/2} \frac{1-2j/m}{\cos(\pi j/m)} \right] \quad \text{if } m \text{ is even,} \qquad (C.32)$$

$$\sigma_m = m \left[ \frac{1}{2} + \sum_{j=1}^{(m-1)/2} \frac{1-2j/m}{\cos(\pi j/m)} \right] \quad \text{if } m \text{ is odd.} \qquad (C.33)$$

As a consequence, we get the variance, $\text{Var}_0\{\ell_t\} = L_2(t|0) - [L_1(t|0)]^2 = 4Dt(\sigma_m - m^2/\pi)$, and thus the prefactor $v_\alpha$ from Eq. (15) reads

$$v_\alpha = 2 \frac{\sigma_m - m^2/\pi}{m(1 - 2/\pi)}. \qquad (C.34)$$

# Appendix D. Mathematical basis for the EFS approach

In this Appendix, we derive the main formulas needed for the implementation of the EFS approach. We consider a sector of angle $\alpha$ and radius $\epsilon$, which is defined in polar coordinates as: $\Omega = \{(r, \theta) : 0 < r < \epsilon, 0 < \theta < \alpha\}$ (see Fig. 2).

## Appendix D.1. Mean escape time

We first look at the mean first-passage time (MFPT) to the absorbing arc $\partial\Omega_D$ of the sector, when the particle started from a point $\mathbf{x}_0 \in \Omega: T(\mathbf{x}_0) = \mathbb{E}_{\mathbf{x}_0}\{\tau\}$. This function satisfies:

$$\Delta T = -\frac{1}{D}, \qquad T_{r=\epsilon} = 0, \qquad (\partial_{\theta} T)_{\theta \in \{0, \alpha\}} = 0, \qquad (D.1)$$

which admits an explicit exact solution:

$$\mathbb{E}_{\mathbf{x}_0}\{\tau\} = T = \frac{\epsilon^2 - r_0^2}{4D}. \qquad (D.2)$$

Indeed, as the segments $\Gamma$ are considered as reflecting here, they have no effect onto radial displacements so that one retrieves the MFPT to the boundary of a disk.

## Appendix D.2. Mean escape position

Next, we need to find the mean escape position on the arc (i.e., the angle $\theta_\tau$). For this purpose, one can compute the Green's function satisfying mixed Dirichlet/Neumann boundary conditions:

$$-\Delta G(\mathbf{x}, \mathbf{x}_0) = \delta(\mathbf{x} - \mathbf{x}_0), \qquad G|_{\partial\Omega_D} = 0, \qquad \partial_n G|_{\Gamma} = 0. \qquad (D.3)$$
As this computation is standard, we just sketch the main steps. We search the Green's function in the form

$$G(\boldsymbol{x}, \boldsymbol{x}_0) = g_0(r, r_0) + \sum_{n=1}^{\infty} \cos(\nu_n \theta) \cos(\nu_n \theta_0) g_n(r, r_0), \qquad (D.4)$$

where $\nu_n = \pi n / \alpha$, and $g_n(r, r_0)$ are unknown radial functions that need to be determined from the equation

$$(\partial_r^2 + r^{-1}\partial_r - \nu_n^2 r^{-2})g_n(r, r_0) = -\frac{2}{\alpha r_0}\delta(r - r_0). \qquad (D.5)$$

For $n > 0$, we search for a solution in a standard way:

$$g_n(r, r_0) = \begin{cases} A_n r^{\nu_n} & (r < r_0), \\ C_n((r/\epsilon)^{\nu_n} - (r/\epsilon)^{-\nu_n}) & (r > r_0), \end{cases} \qquad (D.6)$$

where $A_n$ and $C_n$ are unknown coefficients. This form ensures the regular behavior at $r = 0$ and vanishing of $G$ at $r = \epsilon$. Matching the two parts at $r = r_0$ (by requiring the continuity of $G$ and the drop of its derivative), we get

$$C_n = \frac{a_n r_0}{2\nu_n} (r_0/\epsilon)^{\nu_n}, \qquad A_n = \frac{C_n}{r_0^{\nu_n}} [(r_0/\epsilon)^{\nu_n} - (r_0/\epsilon)^{-\nu_n}], \qquad (D.7)$$

for any $n > 0$. In turn, for $n = 0$, we search $g_0(r, r_0)$ as

$$g_0(r, r_0) = \begin{cases} A_0 & (r < r_0), \\ C_0 \ln(r/\epsilon) & (r > r_0). \end{cases} \qquad (D.8)$$

After matching these solutions at $r = r_0$, we get

$$g_0(r, r_0) = \frac{1}{\alpha} \times \begin{cases} \ln(\epsilon/r_0) & (r < r_0), \\ \ln(\epsilon/r) & (r > r_0). \end{cases} \qquad (D.9)$$

As a consequence, we have

$$G(\boldsymbol{x}, \boldsymbol{x}_0) = g_0(r, r_0) + \sum_{n=1}^{\infty} \frac{1}{\pi n} \cos(\nu_n \theta) \cos(\nu_n \theta_0) \times \begin{cases} [(r/r_0)^{\nu_n} - (rr_0/\epsilon^2)^{\nu_n}] & (r < r_0), \\ [(r_0/r)^{\nu_n} - (rr_0/\epsilon^2)^{\nu_n}] & (r > r_0). \end{cases} \qquad (D.10)$$

Knowing the Green's function, we evaluate the harmonic measure density on the arc:

$$\omega(\theta|\boldsymbol{x}_0) = -(\partial_n G)_{r=R} = \frac{1}{\alpha R} \left( 1 + 2 \sum_{n=1}^{\infty} \cos(\nu_n \theta) \cos(\nu_n \theta_0) (r_0/\epsilon)^{\nu_n} \right). \qquad (D.11)$$

Using the geometric series formula, one has

$$\sum_{n=1}^{\infty} \cos(n\beta)\zeta^n = \zeta \frac{\cos\beta - \zeta}{1 - 2\zeta\cos\beta + \zeta^2}, \qquad (D.12)$$
so that

$$ \omega(\theta|\mathbf{x}_0) = \frac{1}{\alpha R} \left( 1 + \zeta \frac{\cos(\pi(\theta - \theta_0)/\alpha) - \zeta}{1 - 2\zeta \cos(\pi(\theta - \theta_0)/\alpha) + \zeta^2} + \zeta \frac{\cos(\pi(\theta + \theta_0)/\alpha) - \zeta}{1 - 2\zeta \cos(\pi(\theta + \theta_0)/\alpha) + \zeta^2} \right), \quad (D.13) $$

where $\zeta = (r_0/\varepsilon)^{\pi/\alpha}$.

In particular, we can compute the mean escape angle:

$$ \mathbb{E}_{\mathbf{x}_0}\{\theta_{\tau}\} = \varepsilon \int_{0}^{\alpha} d\theta \, \theta \, \omega(\theta|\mathbf{x}_0) = \frac{\alpha}{2} \left( 1 + 4 \sum_{n=1}^{\infty} \frac{(-1)^n - 1}{\pi^2 n^2} \cos(\nu_n \theta_0) (r_0/\varepsilon)^{\nu_n} \right). \quad (D.14) $$

For instance, if $\theta_0 = \alpha/2$, one gets $\mathbb{E}_{\mathbf{x}_0}\{\theta_{\tau}\} = \alpha/2$ as expected.

## Appendix D.3. Mean acquired boundary local time

Finally, we aim at evaluating the mean boundary local time $\ell_\tau$ acquired up to the escape moment $\tau$. Following [72], we consider the joint probability density function of $(\mathbf{X}_\tau, \ell_\tau, \tau)$, denoted as $j_D(\mathbf{x}, \ell, t|\mathbf{x}_0)$, with $\mathbf{x} \in \partial\Omega_D$. On the one hand, the integral of this quantity over $t$ and $\mathbf{x} \in \partial\Omega_D$ yields the (marginal) probability density function of $\ell_\tau$:

$$ \rho_D(\ell|\mathbf{x}_0) = \int_{\partial\Omega_D} d\mathbf{x} \int_0^\infty dt \, j_D(\mathbf{x}, \ell, t|\mathbf{x}_0). \quad (D.15) $$

On the other hand, multiplying this joint PDF by $e^{-q\ell}$ and integrating over $\ell$ yields the probability flux density onto $\partial\Omega_D$ in the presence of a partially reactive boundary $\Gamma$:

$$ j_q(\mathbf{x}, t|\mathbf{x}_0) = \int_0^\infty d\ell \, e^{-q\ell} j_D(\mathbf{x}, \ell, t|\mathbf{x}_0). \quad (D.16) $$

Its integral over $\mathbf{x} \in \partial\Omega_D$ yields the probability density of the FPT to $\partial\Omega_D$ in the presence of a partially reactive boundary $\Gamma$:

$$ J_q(t|\mathbf{x}_0) = \int_{\partial\Omega_D} d\mathbf{x} \int_0^\infty d\ell \, e^{-q\ell} j_D(\mathbf{x}, \ell, t|\mathbf{x}_0). \quad (D.17) $$

In particular, its integral over $t$ is the splitting probability, i.e., the probability of hitting $\partial\Omega_D$ before reacting on $\Gamma$:

$$ \tilde{J}_q(0|\mathbf{x}_0) = \int_0^\infty dt \, J_q(t|\mathbf{x}_0). \quad (D.18) $$

We recall that the splitting probability $\tilde{J}_q(0|\mathbf{x}_0)$ satisfies:

$$ \Delta \tilde{J}_q(0|\mathbf{x}_0) = 0, \quad \tilde{J}_q(0|\mathbf{x}_0) = 1 \quad \text{on } \partial\Omega_D, \quad (\partial_n + q) \tilde{J}_q(0|\mathbf{x}_0) = 0 \quad \text{on } \Gamma. \quad (D.19) $$
Comparing these expressions with Eq. (D.15), we finally get

$$ \tilde{J}_q(0|\mathbf{x}_0) = \int_0^\infty d\ell e^{-q\ell} \rho_D(\ell|\mathbf{x}_0) = \mathbb{E}_{\mathbf{x}_0}\{e^{-q\ell_\tau}\}, \qquad (D.20) $$

i.e., the splitting probability $\tilde{J}_q(0|\mathbf{x}_0)$ is the moment-generating function of $\ell_\tau$. Note that this relation could alternatively be deduced from spectral expansions derived in [72], based on the Steklov-Dirichlet spectral problem. As a consequence, $\tilde{J}_q(0|\mathbf{x}_0)$ determines all positive integer-order moments of $\ell_\tau$. For instance, the mean acquired boundary local time is

$$ \mathbb{E}_{\mathbf{x}_0}\{\ell_\tau\} = -\lim_{q \to 0} \partial_q \tilde{J}_q(0|\mathbf{x}_0). \qquad (D.21) $$

In the limit $q \to 0$, we search $\tilde{J}_q(0|\mathbf{x}_0)$ as a formal expansion:

$$ \tilde{J}_q(0|\mathbf{x}_0) = v_0(\mathbf{x}_0) - qv_1(\mathbf{x}_0) + O(q^2). \qquad (D.22) $$

According to Eq. (D.19), functions $v_0$ and $v_1$ satisfy

$$ \Delta v_0 = 0, \quad v_0|_{\partial\Omega_D} = 1, \quad \partial_n v_0|_{\Gamma} = 0, \qquad (D.23) $$

and

$$ \Delta v_1 = 0, \quad v_1|_{\partial\Omega_D} = 0, \quad \partial_n v_1|_{\Gamma} = v_0. \qquad (D.24) $$

One sees that $v_0 \equiv 1$, whereas $v_1$ can be found as

$$ v_1(\mathbf{x}_0) = \int_{\Gamma} d\mathbf{x} \, G(\mathbf{x}, \mathbf{x}_0) v_0(\mathbf{x}), \qquad (D.25) $$

where $G(\mathbf{x}, \mathbf{x}_0)$ is the Green's function given by Eq. (D.10). Evaluating this integral, we get

$$ v_1(r_0, \theta_0) = 2\frac{\varepsilon - r_0}{\alpha} + \frac{2\varepsilon}{\alpha} \sum_{n=1}^{\infty} \frac{(1 + (-1)^n)}{\nu_n^2 - 1} \cos(\nu_n \theta_0) \left( \frac{r_0}{\varepsilon} - (r_0/\varepsilon)^{\nu_n} \right). \qquad (D.26) $$

We can therefore identify $v_1$ with $\mathbb{E}_{\mathbf{x}_0}\{\ell_\tau\}$. Using the summation formula (see Table 2 from [85]), we can compute explicitly the first sum that yields

$$ \mathbb{E}_{\mathbf{x}_0}\{\ell_\tau\} = \frac{2\varepsilon}{\alpha} - r_0 \frac{\cos(\theta_0) + \cos(\alpha - \theta_0)}{\sin(\alpha)} - \frac{2\varepsilon}{\alpha} \sum_{n=1}^{\infty} \frac{(1 + (-1)^n)}{\nu_n^2 - 1} \cos(\nu_n \theta_0) (r_0/\varepsilon)^{\nu_n}. \qquad (D.27) $$

We can simplify this expression as

$$ \mathbb{E}_{\mathbf{x}_0}\{\ell_\tau\} = \frac{2\varepsilon}{\alpha} - r_0 \frac{\cos(\theta_0) + \cos(\alpha - \theta_0)}{\sin(\alpha)} - \frac{4\varepsilon}{\alpha} \sum_{n=1}^{\infty} \frac{\cos(2\nu_n\theta_0)}{4\nu_n^2 - 1} (r_0/\varepsilon)^{2\nu_n}. \qquad (D.28) $$

For instance, if the starting point $\mathbf{x}_0$ is located on the vertex (i.e., $r_0 = 0$), we get $\mathbb{E}_0\{\ell_\tau\} = 2\varepsilon/\alpha$, i.e., it increases as $\alpha$ decreases.
In the limit $\alpha \to 2\pi$, one needs to treat separately the diverging contributions from the second term and the first term of the sum. Setting $\alpha = 2\pi - \epsilon$ and evaluating the limit $\epsilon \to 0$, we get for $\alpha = 2\pi$:

$$ \mathbb{E}_{x_0}\{\ell_{\tau}\} = \frac{\varepsilon - r_0}{\pi} - r_0 \frac{\cos(\theta_0) \ln(r_0/\varepsilon)}{\pi} + \frac{2\varepsilon}{\pi} \sum_{n=2}^{\infty} \frac{\cos(n\theta_0)}{n^2 - 1} \left(r_0/\varepsilon - (r_0/\varepsilon)^n\right). \quad (D.29) $$

Using the identity

$$ \sum_{n=2}^{\infty} \frac{\cos(n\theta_0)}{n^2 - 1} = \frac{1}{2} + \frac{\cos(\theta_0)}{4} - \sin(\theta_0) \frac{\pi - \theta_0}{2}, \quad (D.30) $$

we get another representation for $\alpha = 2\pi$:

$$ \mathbb{E}_{x_0}\{\ell_{\tau}\} = \frac{\varepsilon}{\pi} + \frac{r_0}{\pi} \left( \frac{\cos(\theta_0)}{2} - (\pi - \theta_0) \sin(\theta_0) - \cos(\theta_0) \ln(r_0/\varepsilon) \right) - \frac{2\varepsilon}{\pi} \sum_{n=2}^{\infty} \frac{\cos(n\theta_0)}{n^2 - 1} (r_0/\varepsilon)^n. \quad (D.31) $$

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

[2] Itô K and McKean H P 1965 *Diffusion Processes and Their Sample Paths* (Berlin: Springer-Verlag)

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

[4] McKean H P Brownian local time, 1975 *Adv. Math.* **15** 91

[5] Borodin A N and Salminen P 1996 *Handbook of Brownian Motion: Facts and Formulae* 2nd ed (Basel-Boston-Berlin: Birkhäuser Verlag)

[6] Majumdar S N Brownian functionals in physics and computer science, 2007 *Curr. Sci.* **89** 2076

[7] Grebenkov D S Paradigm shift in diffusion-mediated surface phenomena, 2020 *Phys. Rev. Lett.* **125** 078102

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

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

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

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

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

[13] Grebenkov D S Depletion of resources by a population of diffusing species, 2022 *Phys. Rev. E* **105** 054402

[14] Grebenkov D S Diffusion-controlled reactions with non-markovian binding/unbinding kinetics, 2023 *J. Chem. Phys.* **158** 214111

[15] Bressloff P C Encounter-based reaction-subdiffusion model i - surface adsorption and the local time propagator, 2023 *J. Phys. A: Math. Theor.* **56** 435004

[16] Bressloff P C Encounter-based reaction-subdiffusion model ii: partially absorbing traps and the occupation time propagator, 2023 *J. Phys. A: Math. Theor.* **56** 435005

[17] Bressloff P C Close encounters of the sticky kind: Brownian motion at absorbing boundaries, 2023 *Phys. Rev. E* **107** 064121

[18] Felici M, Filoche M and Sapoval B Diffusional screening in the human pulmonary acinus, 2003 *J. Appl. Physiol.* **94** 2010

[19] Weibel E R 1984 *The Pathway for Oxygen: Structure and Function in the Mammalian Respiratory System* (Cambridge, Massachusetts and London, England: Harvard University Press)
[20] Sapoval B, Filoche M and Weibel E R Smaller is better – but not too small: A physical scale for the design of the mammalian pulmonary acinus, 2002 *Proc. Natl. Acad. Sci. U.S.A.* **99** 10411

[21] Sapoval B, Kang M Y and Dinh-Xuan A T Modeling of gas exchange in the lungs, 2021 *Compr. Physiol.* **11** 1289

[22] Filoche M, Grebenkov D S, Andrade J S and Sapoval B Passivation of irregular surfaces accessed by diffusion, 2008 *Proc. Natl. Acad. Sci. U.S.A.* **105** 7636

[23] Bond G C 1987 *Heterogeneous Catalysis: Principles and Applications* (Oxford: Clarendon Press)

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

[25] Gheorghiu S and Coppens M O Optimal bimodal pore networks for heterogeneous catalysis, 2004 *AIChE* **50** 812

[26] Sapoval B, Andrade J S and Filoche M Catalytic effectiveness of irregular interfaces and rough pores: the 'land surveyor approximation', 2001 *Chem. Eng. Sci.* **56** 5011

[27] Filoche M, Sapoval B and Andrade J S Deactivation dynamics of rough catalytic surfaces, 2005 *AIChE* **51** 998

[28] de Levie R The influence of surface roughness of solid electrodes on electrochemical measurements, 1965 *Electrochim. Acta* **10** 113

[29] Pajkossy T Electrochemistry at fractal surfaces, 1991 *J. Electroanal. Chem* **300** 1

[30] Sapoval B, Gutfraind R, Meakin P, Keddam M and Takenouti H Equivalent-circuit, scaling, random-walk simulation, and an experimental study of self-similar fractal electrodes and interfaces, 1993 *Phys. Rev. E* **48** 3333

[31] de Levie R Fractals and rough electrodes, 1990 *J. Electroanal. Chem* **281** 1

[32] Halsey T C and Leibig M The double layer impedance at a rough surface: Theoretical results, 1992 *Ann. Phys.* **219** 109

[33] Gabdoulline R R and Wade R C Biomolecular diffusional association, 2002 *Curr. Opin. Struct. Biol.* **12** 204

[34] Basak S, Sengupta S and Chattopadhyay K Understanding biochemical processes in the presence of sub-diffusive behavior of biomolecules in solution and living cells, 2019 *Biophys. Rev.* **11** 851

[35] Lyu Y, An L, Zeng H, Zheng F, Guo J, Zhang P, Yang H and Li H First-passage time analysis of diffusion-controlled reactions in single-molecule detection, 2023 *Talanta* **260** 124569

[36] Collins F C and Kimball G E Diffusion-controlled reaction rates, 1949 *J. Colloid Sci.* **4** 425

[37] Berg H C and Purcell E M Physics of chemoreception, 1977 *Biophys. J.* **20** 193

[38] Sano H and Tachiya M Partially diffusion-controlled recombination, 1979 *J. Chem. Phys.* **71** 1276

[39] Brownstein K R and Tarr C E Importance of classical diffusion in nmr studies of water in biological cells, 1979 *Phys. Rev. A* **19** 2446

[40] Weiss G H Overview of theoretical models for reaction rates, 1986 *J. Stat. Phys.* **42** 3

[41] Bressloff P C, Earnshaw B A and Ward M J Diffusion of protein receptors on a cylindrical dendritic membrane with partially absorbing traps, 2008 *SIAM J. Appl. Math.* **68** 1223

[42] Grebenkov D S Diffusion-controlled reactions: An overview, 2023 *Molecules* **28** 7570

[43] Galanti M, Fanelli D, Angioletti-Uberti S, Ballauff M, Dzubiella J and Piazza F Reaction rate of a composite core-shell nanoreactor with multiple nanocatalysts, 2016 *Phys. Chem. Chem. Phys.* **18** 20758

[44] Sapoval B General formulation of laplacian transfer across irregular surfaces, 1994 *Phys. Rev. Lett.* **73** 3314

[45] Bénichou O and Voituriez R From first-passage times of random walks in confinement to geometry-controlled kinetics, 2014 *Phys. Rep.* **539** 225

[46] Bénichou O, Chevalier C, Klafter J, Meyer B and Voituriez R Geometry-controlled kinetics, 2010 *Nat. Chem.* **2** 472

[47] Grebenkov D S, Metzler R and Oshanin G Strong defocusing of molecular reaction times results from an interplay of geometry and reaction control, 2018 *Commun. Chem.* **1** 96

[48] Varadhan S R S and Williams R J Brownian motion in a wedge with oblique reflection, 1985
Commun. Pure Appl. Math. 38 405

[49] Williams R J Local time and excursions of reflected brownian motion in a wedge, 1987 PRIMS 23
297

[50] DeBlassie R D The first exit time of a two-dimensional symmetric stable process from a wedge, 1990 Ann. Probab. 18 1034

[51] Comtet A and Desbois J Brownian motion in wedges, last passage time and the second arc-sine law, 2003 J. Phys. A: Math. Gen. 36 L255

[52] Lenzi E K, Evangelista L R, Lenzi M K and da Silva L R Fokker-planck equation in a wedge domain: Anomalous diffusion and survival probability, 2009 Phys. Rev. E 80 021131

[53] Bras P and Kohatsu-Higa A Simulation of reflected brownian motion on two dimensional wedges, 2023 Stoch. Proc. Appl. 156 349

[54] Muller M E Some continuous monte carlo methods for the dirichlet problem, 1956 Ann. Math. Stat. 27 569

[55] Zhou Y, Cai W and Hsu E Computation of the local time of reflecting brownian motion and the probabilistic representation of the neumann problem, 2017 Commun. Math. Sci. 15 237

[56] Schumm R D and Bressloff P C A numerical method for solving snapping out brownian motion in 2d bounded domains, 2023 J. Comput. Phys. 493 112479

[57] Ye Y, Chaigneau A and Grebenkov D S Escape-from-a-layer approach for simulating the boundary local time in euclidean domains, 2024 J. Comput. Phys. Accepted; preprint available at arXiv:2411.10220

[58] Akkermans E, Bénichou O, Dunne G V, Teplyaev A and Voituriez R Spatial log-periodic oscillations of first-passage observables in fractals, 2012 Phys. Rev. E 86 061125

[59] Ye Y, Chaigneau A and Grebenkov D S First-passage times to a fractal boundary: local persistence exponent and its log-periodic oscillations, 2025 Phys. Rev. E 111 014153

[60] Grebenkov D S Probability distribution of the boundary local time of reflected brownian motion in euclidean domains, 2019 Phys. Rev. E 100 062110

[61] Dy D L L and Esguerra P First-passage time distribution for diffusion in a planar wedge, 2008 Phys. Rev. E 78 062101

[62] Chupeau M, Bénichou O and Majumdar S Survival probability of a brownian motion in a planar wedge of arbitrary angle, 2015 Phys. Rev. E 91 032106

[63] Le Vot F, Yuste S B, Abad E and Grebenkov D S First-encounter time of two diffusing particles in confinement, 2020 Phys. Rev. E 102 032118

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

[65] Levitin M and Parnovski L On the principal eigenvalue of a robin problem with a large parameter, 2008 Math. Nachr. 281 272

[66] Khalile M and Pankrashkin K Eigenvalues of robin laplacians in infinite sectors, 2018 Math. Nachr. 291 928

[67] Pankrashkin K and Vogel M Asymptotics of robin eigenvalues on sharp infinite cones, 2023 J. Spectr. Theory 13 201

[68] Lyalinov M A A comment on eigenfunctions and eigenvalues of the laplace operator in an angle with robin boundary conditions, 2021 J. Math. Sci. 252 646

[69] Redner S 2001 A Guide to First Passage Processes (Cambridge University Press)

[70] Metzler R, Oshanin G and Redner S (eds) 2014 First-Passage Phenomena and Their Applications (Singapore: World Scientific Press)

[71] Masoliver J 2018 Random Processes: First-passage and Escape (World Scientific)

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

[73] Mandelbrot B 1982 The Fractal Geometry of Nature (San Francisco: Freeman)

[74] Grebenkov D S, Lebedev A A, Filoche M and Sapoval B Multifractal properties of the harmonic measure on koch boundaries in two and three dimensions, 2005 Phys. Rev. E 71 056121

[75] Grebenkov D S Scaling properties of the spread harmonic measures, 2006 Fractals 14 231
[76] Grebenkov D S Analytical representations of the spread harmonic measure density, 2015 *Phys. Rev. E* **91** 052108

[77] Chaigneau A and Grebenkov D S A numerical study of the generalized steklov problem in planar domains, 2024 *J. Phys. A: Math. Theor.* **57** 445201

[78] Traytak S D The diffusive interaction in diffusion-limited reactions: steady-state case, 1992 *Chem. Phys. Lett.* **197** 247

[79] Grebenkov D S What makes a boundary less accessible, 2005 *Phys. Rev. Lett.* **95** 200602

[80] Andrade J S, Araújo A D, Filoche M and Sapoval B Screening effects in flow through rough channels, 2007 *Phys. Rev. Lett.* **98** 194101

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

[82] Bray A J, Majumdar S N and Schehr G Persistence and first-passage properties in non-equilibrium systems, 2013 *Adv. Phys.* **62** 225

[83] Levernier N, Dolgushev M, Bénichou O, Voituriez R and Guérin T Survival probability of stochastic processes beyond persistence exponents, 2019 *Nat. Commun.* **10** 2990

[84] Carslaw H S and Jaeger J C 1959 *Conduction of Heat in Solids* 2nd ed (Oxford University Press)

[85] Grebenkov D S A physicist's guide to explicit summation formulas involving zeros of bessel functions and related spectral sums, 2021 *Rev. Math. Phys.* **33** 2130002