text:
\section{Introduction}\label{sec:intro}

We consider the problem of describing the full arrival time distribution
of diffusing particles to complex absorbing sets in planar regions
$\RR^2\setminus\Omega$ (see Fig.~\ref{Fig:schematic}). For a particle
released from location $\mathbf{x}_0\in\RR^2\setminus\Omega$, the central
quantities of interest are its occupation density $p(\mathbf{x},t;\mathbf{x}_0)$ and
its survival probability
$P(t;\mathbf{x}_0)=\int_{\RR^2\setminus\Omega}p(\mathbf{x},t;\mathbf{x}_0)\,d\mathbf{x}$ together
with the dependence on the number and location of absorbing sites. Our
new contributions are several new efficient numerical and asymptotic
methods to rapidly determine these quantities in the presence of complex
configurations of targets (cf.~Fig.~\ref{Fig:schematic}).

\begin{figure}[htbp]
\centering
\includegraphics[width=0.45\textwidth]{introSchematic.pdf}
\caption{Schematic of planar Brownian motion to a collection of targets $\Omega := \cup_{j=1} \mathcal{A}_j$. The shape of each target can be general and the boundary of each target can be combinations of absorbing and reflecting sections.  \label{Fig:schematic} }
\end{figure}

The general problem takes the form of a diffusion equation in $\mathbf{x} \in
\RR^2\setminus\Omega$ where $\Omega := \cup_{j=1} \mathcal{A}_j$ is a collection of targets. We solve for the probability $p(\mathbf{x},t;\mathbf{x}_0)$ of a particle originating at $\mathbf{x}_0$ being free at position $\mathbf{x}$ at time $t$. This probability density solves the exterior parabolic problem
\begin{subequations}\label{eqn:IntroP}
\begin{gather}
  \label{eqn:IntroP_a}
  \frac{\partial p}{\partial t} = D \Delta p, \qquad
    \mathbf{x} \in \RR^2\setminus\Omega, \quad t>0; \qquad
    p(\mathbf{x},0) = \delta (\mathbf{x}-\mathbf{x}_0), \quad \mathbf{x} \in \RR^2\setminus\Omega;\\[5pt]
  \label{eqn:IntroP_b}
  p = 0, \quad \mathbf{x} \in \partial\Omega_a; \qquad
  D\nabla p \cdot {\mathbf{n}} = 0, \quad \mathbf{x} \in \partial\Omega_r ,
\end{gather}
\end{subequations}$\!$
where $D$ is the diffusivity of the particle and $\Omega$ is a
subset of $\RR^2$. The boundary $\pOm$ is partitioned into an
absorbing set $\Omega_a$ and its impermeable complement $\Omega_r$ where reflecting boundary conditions are applied. We choose ${\mathbf{n}}$, the normal to the surface $\pOm$, to point
into the bulk. Some key quantities of interest for which rapid and
accurate determination is desirable include the fluxes over each target,
the splitting probabilities (likelihood particle encounters a certain
target first) and the arrival time distribution. The survival
probability
\begin{align*}
  P(t;\mathbf{x}_0) = \int_{\RR^2\setminus\Omega} p(\mathbf{x},t;\mathbf{x}_0) \, d\mathbf{x},
\end{align*}
is another important quantity obtained, together with its dependence on the number and location of targets.

First passage time problems and their variants appear in a variety of
disparate applications from cellular biology~\cite{Miles2020,Bressloff2021a,Bressloff2021b,BressNewby2013}, ecology~\cite{Kurella2015, Schenker2015,
Miller2017, Stepien2020}, and electrostatics~\cite{Chapman2015}. A few
comprehensive survey references can be found here~\cite{redner2001guide,
FPPA2014, berg, HolcmanReview2014}. Preceding works on the theory of
first arrival times to small absorbing sites have largely focussed on
determining moments such as the mean first passage time
(MFPT)~\cite{Pillay2010, ChevWard2010, HolcmanReview2014, LIW2021_a,
IW2021_b, IW2021_c, IW2021_d} and in some cases the
variance~\cite{Kurella2015,LWB2017,LTS2016}. In the present scenario of
an unbounded domain, these moments are not finite~\cite{redner2001guide,
FPPA2014} and we must therefore seek the full distribution of arrival
times. In the scenario where diffusion occurs in a bounded domain, the full arrival time distributions to small absorbing targets have been considered in two \cite{LTS2016,Bressloff2021a} and three \cite{Isaacson2013,Bressloff2021b} dimensions.

In the planar (unbounded) scenario considered here, capture is guaranteed; however,
it may occur over very long timescales. This becomes apparent from the arrival time distribution $\chi_0(t)$ for a particle of diffusivity $D=1$, initially at distance $R$ from a target disc of unit radius centered at the origin. This distribution $\chi_0(t)$ and its large time behavior (see \cite{redner_2001} and Appendix \ref{sec:PlanarDist}) are given by
\begin{align}
  \label{eq:dist2Dreinsertion}
  \chi_0(t) = \frac{2}{\pi}\int_{0}^{\infty}\left[
    \frac{J_0(\omega)Y_0(\omega R)-J_0(\omega R)Y_0(\omega)}
    {Y_{0}^{2}(\omega)+J_{0}^{2}(\omega)}\right]\omega
    e^{-\omega^2 t} \, d\omega =
  \mathcal{O}\left(\frac{1}{t\log^2 t}\right), \qquad t\to\infty,
\end{align}
where $J_0(z)$, $Y_0(z)$ are Bessel functions. The slow rate of decay in
the tail of this distribution reveals that very long arrival times are
typical (see Fig.~\ref{fig:2DarrivalCDF_a}). For example, when $R = 10$,
a particle with diffusivity $D=1$ still has an approximately $20\%$
chance of being free after $t = 10^{10}$. Equilibrium quantities (e.g.~splitting probabilities) therefore emerge on timescales that may not necessarily be the most biologically meaningful. For example, in applications such as a moth's search for a mate~\cite{Stepien2020},
or cellular signaling where a downstream event initializes as soon as a
molecule reaches a receptor~\cite{Miles2020}, the
statistics of particles which reach the target first are of most
interest. These \emph{extreme statistics} are governed by the behavior
of $P(t;\mathbf{x}_0)$ for $t\ll1$ \cite{Lawley2020} and so it is necessary to
have methodologies for determining full distributions of arrival time
statistics.

In the present work, we outline several methods to solve
\eqref{eqn:IntroP}. First, in Sec.~\ref{sec:asy} we apply a Laplace
transform $\hat{p}(\mathbf{x},s)= \int_{t=0}^\infty p(\mathbf{x},t) e^{-st}\,dt$ to
\eqref{eqn:IntroP} to arrive at an elliptic problem of modified
Helmholtz type
\begin{subequations}\label{eq:LT2D}
\begin{gather}
  \label{eq:LT2D_b}
  D\Delta\hat{p} - s\, \hat{p} = -\delta(\mathbf{x}-\mathbf{x}_0),\qquad
    \mathbf{x}\in\RR^2\setminus\Omega;\\[5pt]
  \label{eq:LT2D_a}
  \hat{p} = 0 \quad \mathbf{x} \in \Omega_a; \qquad
    \nabla \hat{p} \cdot \textbf{n} = 0 \quad \mathbf{x} \in \Omega_r.
\end{gather}
\end{subequations}$\!$
In the limit of well-separated traps, we solve the resulting transform
problem in terms of an asymptotic expansion where solutions are obtained
in terms of modified Helmholtz Green's function. This methodology was
originally developed in~\cite{LTS2016} and recently applied
in~\cite{Bressloff2021a,Bressloff2021b} to the study of first passage
times of particles with resetting. The Laplace transform is inverted
numerically by Talbot quadrature~\cite{Abate2006} resulting in a hybrid
numerical-asymptotic method. In Sec.~\ref{sec:BIE}, we take a similar
approach, but replace the asymptotic solution of~\eqref{eq:LT2D} with a
layer potential representation. This results in a boundary integral
equation that is solved numerically using a collocation method.

In Sec.~\ref{sec:KMC} we develop a particle based kinetic Monte-Carlo
(KMC) method that evaluates solutions of \eqref{eqn:IntroP} by dividing
the sojourn of particles into projection steps where exact solutions are
available \cite{LBS2018,OKMascagni2004}. This offers a rapid, accurate
and easy to implement method for the solution of \eqref{eqn:IntroP} in very
general geometries. In Sec.~\ref{sec:results} we demonstrate the
applicability of these methods on a variety of examples. In particular,
we provide numerical validation of previously derived homogenization
theories and find them to be highly effective in reproducing the arrival
time distributions. We also investigate the time dependent fluxes into the
targets which are observed to converge very slowly to the static splitting
probabilities that describe the relative flux into each target. This suggests that
a relevant physical or biological time scale should be considered before using
 receptor arrival information to make inferences on environmental conditions.

\section{Asymptotic description of arrival times of particles diffusing
in $\RR^2$}
\label{sec:asy}
In this section we use matched asymptotic expansions to derive an
approximation for the density of a particle diffusing in $\RR^2\setminus\Omega$
in the presence of well-separated target sites. The assumption is that
there are $N$ targets $\mathcal{A}_j$ with centers $\{\mathbf{x}_j\}_{j=1}^N$ so that
the collection of target sites are described by
\begin{align}\label{eq:2DAbsSet}
\Omega = \bigcup_{j=1}^N  \mathcal{A}_j \left (\frac{\mathbf{x} -\mathbf{x}_j} {{\displaystyle \varepsilon}} \right ), \qquad \mathbf{x} = (x,y),
\end{align}
where ${\displaystyle \varepsilon}$ is a parameter controlling the extent of the targets and enforces the well-separated condition as ${\displaystyle \varepsilon}\to0$.
The geometry of individual targets can be quite general.

The aim is to solve the solution of the initial-boundary value problem
\eqref{eqn:IntroP} and determine the free probability
$P(t)=\int_{\Omega} p(\mathbf{x},t)\,d\mathbf{x}$ together with the capture time
density $C(t)= -P'(t)$. The first step \cite{LTS2016} in the analysis is to define, for
$s\in\CC$, the Laplace transform $\hat{p}(\mathbf{x},s)= \int_{t=0}^\infty
p(\mathbf{x},t) e^{-st}\,dt$ that solves \eqref{eq:LT2D}. The mixed boundary
conditions \eqref{eq:LT2D_a} indicate that the target boundary may have
a combination of absorbing or reflecting components so that
$\partial\Omega = \Omega_a\cup\Omega_r$. In the absence of the target
set $\Omega$, the solution $\hat{p}$ of \eqref{eq:LT2D} is defined in
terms of the free space modified Helmholtz Green's function
$G_h(\mathbf{x};\boldsymbol{\xi},s)$
\begin{subequations}\label{eqG}
\begin{gather}
\label{eqG_a} D\Delta G_h-s\, G_h=-\delta(\mathbf{x}-\boldsymbol{\xi}), \qquad
  \mathbf{x}\in\RR^2\setminus\{\boldsymbol{\xi}\};\\[5pt]
\label{eqG_b} G_h(\mathbf{x};\boldsymbol{\xi},s)=\frac{1}{2\pi D} K_0(\sqrt{s/D}|\mathbf{x}-\boldsymbol{\xi}|),
  \qquad \int_{\RR^2} G_h(\mathbf{x};\boldsymbol{\xi},s)\,d\mathbf{x}=\frac{1}{s}; \\[5pt]
\label{eqG_c} G_h(\mathbf{x};\boldsymbol{\xi},s) \sim -\frac{1}{2\pi D} \log |\mathbf{x}-\boldsymbol{\xi}| + R_h (s) +  \mbox{{\scriptsize $\mathcal{O}$}}(1), \qquad \mathbf{x}\to\boldsymbol{\xi}.
\end{gather}
\end{subequations}$\!$
Here $R_h(s)$ is the regular part of $G_h(\mathbf{x};\boldsymbol{\xi},s)$ at the source. The small argument asymptotics $K_0(z)\sim-\log(z)+\log2-\gamma_e$ as $z\to0$ give this self-interaction term to be
\begin{equation}
 R_h(s)=\frac{1}{2\pi D}\left(\log2-\gamma_e-\log\sqrt{s/D}\right),
\end{equation}
where $\gamma_e\approx 0.5772$ is the Euler-Mascheroni constant.

In the limit of well-separated absorbers ${\displaystyle \varepsilon}\to0$, we employ a matched
asymptotic analysis to replace each target \eqref{eq:LT2D_a} by
effective singularity conditions. To establish this singularity
condition, the following change of variables is introduced near the
$j^{\mathrm{th}}$ absorber
\begin{equation}
 \mathbf{y}=\frac{\mathbf{x}-\mathbf{x}_j}{{\displaystyle \varepsilon}}, \qquad v(\mathbf{y}) = \hat{p}(\mathbf{x}_j + {\displaystyle \varepsilon} \mathbf{y}).
\end{equation}
In these coordinates, the transformed equation \eqref{eq:LT2D_b} is
$D\Delta_{\mathbf{y}}v - s {\displaystyle \varepsilon}^2 v= - {\displaystyle \varepsilon}^2 \delta(\mathbf{x}-\mathbf{x}_0)$. In addition
to the limit ${\displaystyle \varepsilon}\to0$, we additionally consider the case $s{\displaystyle \varepsilon}^2
\ll1$ which is valid provided $s$ is not too large. The limit
$s\to\infty$ corresponds to $t\to0$ and therefore we cannot expect good
agreement for arbitrarily short times. With these points in mind, we
continue by considering the local solution $v(\mathbf{y}) = v_j (\mathbf{y})+
\mathcal{O}({\displaystyle \varepsilon})$ near the target where $v_j(\mathbf{y})$ satisfies the exterior
problem
\begin{subequations}\label{eqn:LogCap}
\begin{gather}
\Delta_{\mathbf{y}}v_j=0,\qquad \mathbf{y}\in\RR^2/{\mathcal{A}}_j,\\[5pt]
v_j(\mathbf{y})\sim \log|\mathbf{y}|-\log d_j + \frac{ {\textbf p} \cdot \mathbf{y} }{|\mathbf{y}|^2} + \cdots ,\quad |\mathbf{y}|\to\infty.
\end{gather}
\end{subequations}$\!$
Here the parameter $d_j$ is the logarithmic capacitance and depends on
the shape of ${\mathcal{A}}_j$ and the boundary conditions applied to it. In
section \ref{sec:LogCap} we give an overview of many scenarios in which
$d_j$ can be calculated. The behavior of $v_j(\mathbf{y})$ at infinity gives
the matching condition for $\hat{p}(\mathbf{x})$ as $\mathbf{x}$ approaches $\mathbf{x}_j$.
That is,
\begin{align}\label{eq:guage}
\hat{p}(\mathbf{x},s)\sim S_j \nu_j
  \left(\log\left|\frac{\mathbf{x}-\mathbf{x}_j}{{\displaystyle \varepsilon}}\right|-\log d_j\right)=S_j
  \nu_j \log|\mathbf{x}-\mathbf{x}_j|+S_j,\quad \mathbf{x}\to\mathbf{x}_j; \qquad \nu_j = -\frac{1}{\log {\displaystyle \varepsilon} d_j},
\end{align}
where $S_j$ is a strength term to be determined in terms of $(N+1)$ parameters $(s,\bs{\nu})=(s,\nu_1,\ldots,\nu_N)$. Therefore, in the outer region away from targets, we pose the asymptotic expansion
\begin{align*}
  \hat{p}(\mathbf{x};s) = \hat{p}_0(\mathbf{x};s,\bs{\nu}) + {\displaystyle \varepsilon} \,\hat{p}_1(\mathbf{x};s,\bs{\nu}) +  \mbox{{\scriptsize $\mathcal{O}$}}({\displaystyle \varepsilon}), \qquad {\displaystyle \varepsilon}\to0.
\end{align*}
The leading order solution $\hat{p}_0(\mathbf{x};s,\bs{\nu})$ \lq\lq sums-the-logs\rq\rq\ and is accurate to all logarithmic orders. The correction term $\hat{p}_1$, which we do not explicitly determine, describes how target orientation influences capture and can be found following methods outlined in \cite{LTK2015}. The leading order problem satisfies
\begin{subequations}\label{equ2}
\begin{gather}
D\Delta\hat{p}_0 -s\, \hat{p}_0=-\delta(\mathbf{x}-\mathbf{x}_0),\quad \mathbf{x}\in\RR^2/\{\mathbf{x}_1,\ldots,\mathbf{x}_N\},\\[5pt]
\hat{p}_0\sim S_j \nu_j \log|\mathbf{x}-\mathbf{x}_j|+S_j,\quad \mathbf{x}\to\mathbf{x}_j, \qquad j = 1,\ldots, N.
\end{gather}
\end{subequations}$\!$
The solution $\hat{p}_0(\mathbf{x}, s)$ of \eqref{equ2} is described in terms of the modified Helmholtz Green's function \eqref{eqG} as
\begin{equation}
\hat{p}_0(\mathbf{x},s)= G_h(\mathbf{x};\mathbf{x}_0,s)- 2\pi D\sum_{j=1}^N S_j \nu_j G_h(\mathbf{x};\mathbf{x}_j,s).
\label{usol}
\end{equation}
The coefficients $S_j$ can be determined by equating the regular parts in (\ref{equ2}b) and in (\ref{usol}). That is,
\begin{equation}
S_j=
  G_h(\mathbf{x}_j;\mathbf{x}_0,s)- 2\pi D\left[S_j\nu_jR_h(s)+\sum_{\substack{i=1\\
  i\neq j}}^NS_i\nu_iG_h(\mathbf{x}_j;\mathbf{x}_i,s)\right], \qquad j = 1,\ldots,N.
\label{eqA}
\end{equation}
\begin{subequations}\label{sysmain}
In summary, we have that the transform equation \eqref{eq:LT2D} has asymptotic solution $p(\mathbf{x};s) \sim p_0(\mathbf{x};s)+\cdots$ as ${\displaystyle \varepsilon}\to0$ where
$p_0(\mathbf{x};s)$ satisfies
\begin{equation}
\hat{p}_0(\mathbf{x},s)= G_h(\mathbf{x};\mathbf{x}_0,s)- 2\pi D\sum_{j=1}^N S_j \nu_j G_h(\mathbf{x};\mathbf{x}_j,s).
\end{equation}
The strengths $\{S_j\}_{j=1}^N$ satisfy \eqref{eqA} which can be represented in compact matrix form as
\begin{equation}
 (\mathcal{I} + 2\pi D\,  \mathcal{G}_h \mathcal{V})\, \textbf{S}=\mathbf{g}_0, \qquad [\mathcal{G}_h]_{i,j}= \left\{ \begin{array}{rl} R_h(s),\quad i= j, \\[5pt]  G_h(\mathbf{x}_i;\mathbf{x}_j,s),\quad i\neq j. \end{array} \right. \qquad [\mathcal{V}]_{i,j}= \left\{ \begin{array}{rl} \nu_i \quad i= j, \\[5pt]  0,\quad i\neq j, \end{array} \right.
\label{matrixG}
\end{equation}
where $\mathcal{I}\in\RR^{N\times N}$ is the identity,
$\textbf{S}\in\RR^N$ and $ \mathbf{g}_0\in\RR^N$ are given by
\begin{align}
\textbf{S}= [S_1,S_2,\ldots,S_N]^T,\qquad
 \mathbf{g}_0=[G_h(\mathbf{x}_1;\mathbf{x}_0,s), G_h(\mathbf{x}_2;\mathbf{x}_0,s), \ldots,G_h(\mathbf{x}_N;\mathbf{x}_0,s)]^T.
\end{align}
\end{subequations}$\!$
The matrix $\mathcal{G}$ describes the interactions between targets and
their competition for flux while the vector $\mathbf{g}_0$ reflects the
influence of the initial location on each of the targets.  The vector
$\widehat{\mathcal{J}} = 2\pi D \mathcal{V}\textbf{S}(s)$ describes the transformed fluxes through each of the targets and is obtained from solving the linear system \eqref{sysmain}.

At this stage we calculate additional quantities of interest, namely the survival probability $P(t)$ and arrival time distribution $C(t)$. Using equation \eqref{usol} and \eqref{eqG_c}, the Laplace transform $\hat{P}(s)$ of the free probability is given by
\begin{subequations}\label{2DResults}
\begin{equation}\label{2DResults_P}
 \hat{P}(s)=\int_{\Omega}
  \hat{p}_0(\mathbf{x},s)\,d\mathbf{x}= \int_{\Omega} G_h(\mathbf{x};\mathbf{x}_0,s)d\mathbf{x} - 2\pi D\sum_{j=1}^N S_j \nu_j \int_{\Omega} G_h(\mathbf{x};\mathbf{x}_j,s)d\mathbf{x} =\frac{1}{s}\left[1-2\pi D\sum_{j=1}^N \nu_j S_j(s)\right].
\end{equation}
The relationship $C(t)=-P'(t)$, yields that the Laplace transform of the arrival time distribution $C(t)$ is
\begin{equation}\label{2DResults_C}
\hat{C}(s)=-[s\hat{P}(s)-P(0)]=-s\hat{P}(s)+1=2\pi D\sum_{j=1}^N \nu_j S_j(s).
\end{equation}
\end{subequations}$\!$
\begin{figure}[htbp]
\centering
\includegraphics[width = 0.5\textwidth]{TalbotSchematic.pdf}
\caption{A schematic of the Talbot curve \eqref{eq:Talbot} for parameter values $\sigma=0$, $\mu = 2$, $\nu = 0.5$. The red line indicates the singularities along the negative real line arising from the $\sqrt{s}$ singularity. \label{fig:Talbot} }
\end{figure}

\subsection{Inverse Laplace Transform}\label{sec:invLT}

To obtain $P(t)$ and $C(t)$ defined by \eqref{2DResults}, the inverse Laplace transform
\begin{equation}\label{eqn:numInvLT}
P(t) = \frac{1}{2\pi i} \int_{\Gamma_B} e^{st} \hat{P}(s)  \, \textrm{d}s,
\end{equation}
must be evaluated where $\Gamma_B$ is the Bromwich contour $\Gamma_B = \{ \gamma + i y \, | \, -\infty<y<\infty \}$. The parameter $\gamma$ is chosen so that all singularities of $\hat{P}(s)$ lie to the left of Re($s)=\gamma$. In the present scenario associated with diffusive motion, the singularities of $\hat{P}(s)$ lie along the negative real axis due to the branch cut of $\sqrt{s}$. Rapid and effective numerical evaluation of \eqref{eqn:numInvLT} can be achieved by deforming the contour around Re($s)=0$ since the integrand of \eqref{eqn:numInvLT} decays very rapidly for Re($s)<0$. The Talbot contour $\Gamma_T$ is a family of deformations (see Fig.~\ref{fig:Talbot})  to $\Gamma_B$ where
\begin{equation}\label{eq:Talbot}
\Gamma_T = \{ \sigma + \mu (\theta \cot \theta + \nu i \, \theta)  \, | \, -\pi<\theta<\pi \},
\end{equation}
and $\sigma, \mu$ and $\nu$ are parameters that control the curve shape
\cite{Weideman2006,Trefethen2006}. Rapid and accurate evaluation of the
inverse Laplace transform is then achieved by applying the midpoint rule
on this curve.

\subsection{Logarithmic capacitance for various shapes}\label{sec:LogCap}

The asymptotic solution \eqref{2DResults} encodes the geometry of each target $\mathcal{A}_j$ into the logarithmic capacitance $d_j$, determined by the solution of \eqref{eqn:LogCap}. Here we discuss the determination of $d_j$ and briefly recap known results for regular geometries and simple boundary conditions.

\paragraph{Regular geometric shapes} For simple shapes such as circles, ellipses, triangles and squares with all absorbing perimeters, the logarithmic capacitance is known exactly. A list of these quantities, reproduced from \cite{Kurella2015}, is included in Table~\ref{tab:LogCap}.

\begin{table}[H]
    \centering
    \begin{tabular}{c|c}
    \toprule
         {Shape of ${\mathcal{A}}_j$} & Logarithmic capacitance $d_j$  \\
         \midrule
         \midrule
         circle of radius $a$ & $d_j = a$\\[5pt]
         ellipse, semi-axes $a,b$ & $d_j = \frac{a+b}{2}$\\[5pt]
         equilateral triangle, side-length $h$ & $d_j = \frac{\sqrt{3}\Gamma(\frac13)^3 h}{8\pi^2} \approx 0.422 h$\\[5pt]
         isosceles right triangle, side-length $h$ & $d_j = \frac{3^{3/4}\Gamma(1/4)^2 h}{2^{7/2}\pi^{3/2}} \approx 0.476 h$\\[5pt]
         square, side-length $h$ & $d_j = \frac{\Gamma(1/2)^2 h}{4\pi^{3/2}} \approx 0.590 h$\\
         \bottomrule
    \end{tabular}
    \bigskip
    \caption{The logarithmic capacitance of some simple geometries with
    absorbing boundary conditions, reproduced from \cite{Kurella2015}.
    \label{tab:LogCap}}
\end{table}

\paragraph{Partially absorbing disk: Single window} For a circular trap that is absorbing except for the reflecting portion $\theta \in(-\sigma,\sigma)$, the problem \eqref{eqn:LogCap} may be expressed in polar coordinates $(r,\theta)$ as
\begin{subequations}\label{eqn:vc_partial}
\begin{gather}
    \Delta v = 0, \quad r \geq 1, \quad \theta\in(-\pi,\pi);\\[5pt]
    \partial_r v = 0, \qquad r = 1, \quad \theta \in(-\sigma,\sigma); \qquad
    v = 0,\qquad r = 1, \quad \theta \in(\sigma,\pi)\cup (-\pi,-\sigma);\\[5pt]
    v \sim \log r - \log d_c + \mathcal{O}\big(r^{-1}\big), \quad r \to \infty.
\end{gather}
\end{subequations}$\!$
The separable solution of \eqref{eqn:vc_partial} takes the form of the cosine series
\begin{align*}
  v(r,\theta) = \log r + \frac{a_0}{2} +
    \sum_{n=1}^{\infty}\frac{a_n}{r^n}\cos n\theta,
\end{align*}
where the coefficients $\{a_n\}_{n=0}^{\infty}$ satisfy the dual
trigonometric series
\begin{subequations}\label{eqn:dualTrig}
\begin{align}
    \frac{a_0}{2} + \sum_{n=1}^{\infty} a_n \cos n \theta &= 0, \quad \theta\in(\sigma,\pi),\\
    \sum_{n=1}^{\infty} n a_n \cos n \theta &= 1, \quad \theta\in(0,\sigma).
\end{align}
\end{subequations}$\!$
The solution of \eqref{eqn:dualTrig} was determined in \cite{LTK2015} from an integral equation theory which reveals that
\begin{equation}\label{eqn:a0split}
   -\log d_c = \frac{a_0}{2} = \frac{\sqrt{2}}{\pi} \int_0^{\sigma} \frac{u \sin \frac{u}{2}}{\sqrt{\cos u - \cos \sigma}}\, \mathrm{d}u.
\end{equation}
For the half absorbing case ($\sigma = \pi/2$), $a_0 = 2\log2$, while in the singular mostly reflecting limit $\sigma\ll1$, it can be determined that $a_0 \sim -4\log\frac{\sigma}{2}$. In other scenarios, the integral \eqref{eqn:a0split} is readily evaluated by quadrature.

\paragraph{Partially absorbing disk: Multiple windows}
In the scenario of a circular trap with $N$ small absorbing windows of length $\sigma$ centered at points $\{\mathbf{y}_j \}_{j=1}^N$, it was shown in \cite{LTK2015} that as $\sigma\to0$
\begin{subequations}\label{eq:multi_windows}
\begin{equation}\label{eq:multi_windows_a}
-\log d_c = \frac{a_0}{2} \sim -\frac{2}{N} \log\frac{\sigma}{4} - \frac{2}{N^2}\sum_{i=1}^N \log \prod_{\substack{j=1 \\ j\neq i}}^N |\mathbf{y}_i - \mathbf{y}_j|.
\end{equation}
For windows centered at roots of unity $\mathbf{y}_j = \left(\cos \frac{2\pi
j}{N},\sin \frac{2\pi j}{N}\right)$, and for $N\sigma <1$,
\eqref{eq:multi_windows_a} reduces to
\begin{equation}
  \label{eq:multi_windows_b}
  -\log d_c = \frac{a_0}{2} = -\frac{2}{N} \log\frac{\sigma}{4} -
    \frac{2}{N^2}(N\log N)= -\frac{2}{N} \log \frac{\sigma N}{4}.
\end{equation}
\end{subequations}$\!$
\paragraph{Partially absorbing disk: Homogenization limit}
The results \eqref{eq:multi_windows} can be used to identify an
homogenization limit as $N\to \infty$ and $\sigma\to0$. The  absorbing
fraction $f$ is defined through $N\sigma = 2\pi f$ and the homogenized
logarithmic capacitance problem satisfies
\begin{subequations}\label{eqn:homogenization}
\begin{gather}
\label{eqn:homogenization_a} \Delta v_h=0,\quad r>1; \qquad \sigma \frac{\partial v_h}{\partial r} + \kappa(f) v_h = 0, \quad r = 1;\\[5pt]
\label{eqn:homogenization_b} v_h\sim \log r-\log d_c + \cdots ,\quad r\to\infty.
\end{gather}
In the dilute limit $f\ll1$, the homogenized parameters were identified in \cite{LTK2015} to be
\begin{equation}\label{eqn:homogenization_c}
  \kappa(f) = -\frac{\pi f}{\log \frac{\pi f}{2}}, \qquad  \log d_c = \frac{\sigma}{\kappa(f)}\, .
\end{equation}
\end{subequations}$\!$
In Section~\ref{sec:resultsHomogenization}, we show numerical results
that validate this homogenized formulation and demonstrate that it is highly accurate in predicting the arrival time statistics of the full problem.

\paragraph{The logarithmic capacitance for a two trap cluster}
In the case of two circular traps separated by distance $\ell$, it was
derived in \cite{Kurella2015} (see also
\cite{LIW2021_a,IW2021_b,Strieder08,Strieder12}) from an expansion in
bi-polar coordinates that
\begin{equation}\label{eq:TwoTrapCluster}
\log d_c = \frac{1}{2}\log(\ell^2 - 4) - \frac{\beta}{2} + \sum_{k=1}^{\infty}\frac{e^{-k\beta}}{k\cosh(k\beta)}, \qquad \beta = \cosh^{-1} \left(\frac{\ell}{2}\right).
\end{equation}

\paragraph{Numerical evaluation of the logarithmic capacitance for general configurations}

For very general configurations of clusters, the logarithmic capacitance problem \eqref{eqn:LogCap} can be obtained numerically by a boundary integral approach \cite{Dijkstra2008}. Another approach developed in \cite{LIW2021_a,IW2021_b}, and based on works \cite{Chapman2015,trefethen_2018}, is to develop a series solution of \eqref{eqn:LogCap} followed by a least squared method to obtain the unknown coefficients. For the case of $m$ circular absorbers with centers $\{c_j\}_{j=1}^m$, the series takes the form
\begin{equation}\label{eqn:complexSeries}
    v(z) = - \log d_c + \sum_{j=1}^m e_j \log|z - c_j| + \sum_{j=1}^m
    \sum_{k = 1}^n \left( a_{j k} \text{Re}(z-c_j)^{-k} + b_{j k}
    \text{Im} ( z- c_j) ^{-k} \right), \quad z\in\CC.
\end{equation}
The constants $d_c, e_j,a_{jk}, b_{jk}$ are to be determined while $\sum_{j=1}^{m}e_j=1$ enforces the far field behavior $v \sim \log |z| $ as $|z|\to\infty$. A system for the unknown constants is formed by evaluating \eqref{eqn:complexSeries} at a collection of boundary points along which $v = 0$. A similar methodology was used in \cite{LTK2015} to solve a truncated version of the dual trigonometric series \eqref{eqn:dualTrig} by evaluation at a set of boundary values followed by least squared solution.

\section{Boundary integral equation description of arrival times of
particles diffusing in $\RR^2$}
\label{sec:BIE}
An alternative approach to matched asymptotics for
solving~\eqref{eqn:IntroP} is to use an integral equation approach.
Integral equations are a natural choice for unbounded complex domains
such as the one in Fig.~\ref{Fig:schematic} since they easily resolve
complex geometries while automatically satisfying the far-field boundary
conditions. Others have applied integral equation methods to
solve~\eqref{eqn:IntroP} using the full space-time heat
kernel~\cite{gre-li2009} or by discretizing in time and solving the
resulting elliptic PDE with an integral equation
formulation~\cite{kro-qua2010, cha-kre1997}. However, these approaches
have several challenges that include maintaining long time
histories and computing volume integrals. We take a new approach by
solving for the Laplace transformed variable $\hat{p}({\mathbf{x}},s)$ that
satisfies~\eqref{eq:LT2D}. We only consider the case where ${\partial}\Omega$
is absorbing so that the boundary condition is Dirichlet and
homogeneous.

We begin by writing $\hat{p}$ as the sum of a particular and homogeneous
solution of~\eqref{eq:LT2D}
\begin{align*}
  \hat{p}({\mathbf{x}},s) = G_h({\mathbf{x}};{\mathbf{x}}_0,s) + \hat{p}^{H}({\mathbf{x}},s),
\end{align*}
where $G_h$ is the free space modified Helmholtz Green's
function~\eqref{eqG}. Using the boundary condition~\eqref{eq:LT2D_a},
$\hat{p}^H$ satisfies the homogeneous PDE
\begin{subequations}
\label{eqn:homoPDE_BC}
\begin{alignat}{3}
  D \Delta \hat{p}^{H} - s \hat{p}^H &= 0, &&{\mathbf{x}} \in \Omega,
  \label{eqn:homoPDE} \\[5pt]
  \hat{p}^{H}({\mathbf{x}}) &= f({\mathbf{x}}), \quad && {\mathbf{x}} \in {\partial}\Omega,
  \label{eqn:homoBC}
\end{alignat}
\end{subequations}
where $f({\mathbf{x}}) = -G_h({\mathbf{x}};{\mathbf{x}}_0,s)$. We represent the solution
of~\eqref{eqn:homoPDE_BC} with the double-layer potential
\begin{align*}
  \hat{p}^{H}({\mathbf{x}}) = {\mathcal{D}}[\sigma]({\mathbf{x}}) =
    \int_{{\partial}\Omega} \pderiv{}{{{\mathbf{n}}}_{\mathbf{y}}}
    G_{h}({\mathbf{x}};{\mathbf{y}},s) \sigma({\mathbf{y}}) \, ds_{\mathbf{y}}, \quad {\mathbf{x}} \in \Omega,
\end{align*}
where $\sigma$ is an unknown density function. We remind the reader that
the unit normal ${\mathbf{n}}_{\mathbf{y}}$ points into the bulk. To satisfy the boundary
condition~\eqref{eqn:homoBC}, the density function must solve
\begin{align}
  \label{eqn:SKIE}
  f({\mathbf{x}}) = \frac{1}{2} \sigma({\mathbf{x}}) + {\mathcal{D}}[\sigma]({\mathbf{x}}),
    \quad {\mathbf{x}} \in {\partial}\Omega.
\end{align}
A numerical solution of the second-kind integral
equation~\eqref{eqn:SKIE} is formed by discretizing ${\partial}\Omega$ at $N$
quadrature points and approximating the integrals with the trapezoid
rule. The resulting linear system is
\begin{align}
  \label{eqn:discretizedSKIE}
  f({\mathbf{x}}_i) = \frac{1}{2} \sigma({\mathbf{x}}_i) +
    \frac{1}{2\pi} \sum_{j=1}^{N} G_{h}({\mathbf{x}}_i;{\mathbf{x}}_j,s)
    \sigma({\mathbf{x}}_j) \, \Delta s_j, \quad i=1,\ldots,N,
\end{align}
and the diagonal term of this linear system is replaced with the
limiting value
\begin{align*}
  \lim_{\substack{{\mathbf{y}} \rightarrow {\mathbf{x}} \\ {\mathbf{x}} \in {\partial}\Omega}}
    G_h({\mathbf{x}};{\mathbf{y}},s) = -\frac{1}{2}\kappa({\mathbf{x}}),
\end{align*}
where $\kappa({\mathbf{x}})$ is the curvature of ${\partial}\Omega$. The linear
system~\eqref{eqn:discretizedSKIE} is solved with the generalized
minimal residual method (GMRES), and since it is the discretization of a
second-kind integral equation, the number of required iterations is
mesh-independent. This method to solve for $\hat{p}({\mathbf{x}},s)$ is coupled
with the inverse Laplace transform~\eqref{eqn:numInvLT}, where we use
the same Talbot contour illustrated in Fig.~\ref{fig:Talbot}.

The flux at point ${\mathbf{x}} \in {\partial}\Omega$ and time $t$ is ${\mathcal{J}}({\mathbf{x}},t) =
\pderiv{}{{\mathbf{n}}}p({\mathbf{x}},t)$. Since we write $\hat{p}({\mathbf{x}},s)$ as the sum of a
fundamental solution and a homogeneous solution, we compute the flux of
these terms individually, and the flux due to the fundamental solution
is computed analytically. The flux due to the homogeneous solution
$\hat{p}^H$ is
\begin{align}
  \label{eqn:DLPflux}
  \widehat{{\mathcal{J}}}^H({\mathbf{x}};s) = \pderiv{}{{\mathbf{n}}_{\mathbf{x}}} {\mathcal{D}}[\sigma]({\mathbf{x}})
  = \int_{{\partial}\Omega} \pderiv{}{{\mathbf{n}}_{\mathbf{x}}}\pderiv{}{{\mathbf{n}}_{\mathbf{y}}}
  G_h({\mathbf{x}};{\mathbf{y}},s) \sigma({\mathbf{y}}) \, ds_{\mathbf{y}},
\end{align}
which needs to be estimated with quadrature. The trapezoid rule or
odd-even integration~\cite{sid-isr1988} cannot be used to
approximate~\eqref{eqn:DLPflux} because
\begin{align*}
  \pderiv{}{{\mathbf{n}}_{\mathbf{x}}}\pderiv{}{{\mathbf{n}}_{\mathbf{y}}} G_h({\mathbf{x}};{\mathbf{y}},s) \sim
    \|{\mathbf{x}} - {\mathbf{y}}\|^{-2}.
\end{align*}
To formulate the normal derivative of the double-layer potential with a
tractable integrand, we first add and subtract the leading order
asymptotics of $G_h$ described in~\eqref{eqG_c}. That is,
$\widehat{{\mathcal{J}}}^H({\mathbf{x}};s) = I_1 - I_2$ where
\begin{align*}
  I_1 &= \int_{{\partial}\Omega}
    \pderiv{}{{\mathbf{n}}_{\mathbf{x}}}\pderiv{}{{\mathbf{n}}_{\mathbf{y}}}
    \left(G_h({\mathbf{x}};{\mathbf{y}},s) + \frac{1}{2\pi} \log|{\mathbf{x}} - {\mathbf{y}}|\right)
    \sigma({\mathbf{y}}) \, ds_{\mathbf{y}},  \\
  I_2&=\int_{{\partial}\Omega} \pderiv{}{{\mathbf{n}}_{\mathbf{x}}}\pderiv{}{{\mathbf{n}}_{\mathbf{y}}}
    \left(\frac{1}{2\pi} \log|{\mathbf{x}} - {\mathbf{y}}|\right) \sigma({\mathbf{y}}) \, ds_{\mathbf{y}}.
\end{align*}
The singularity of the integrand in $I_1$ behaves as $\|{\mathbf{x}} -
{\mathbf{y}}\|^{-1}$, and odd-even integration can be applied. The integral
$I_2$ is further decomposed as
\begin{align*}
  I_2 &= \frac{1}{2\pi}\int_{{\partial}\Omega}
    \pderiv{}{{\mathbf{n}}_{\mathbf{x}}}\pderiv{}{{\mathbf{n}}_{\mathbf{y}}}\log|{\mathbf{x}} - {\mathbf{y}}|
    (\sigma({\mathbf{y}}) - \sigma({\mathbf{x}})) \, ds_{\mathbf{y}}
    - \frac{\sigma({\mathbf{x}})}{2\pi} \int_{{\partial}\Omega}
      \pderiv{}{{\mathbf{n}}_{\mathbf{x}}}\pderiv{}{{\mathbf{n}}_{\mathbf{y}}} \log|{\mathbf{x}} - {\mathbf{y}}| \, ds_{\mathbf{y}}.
\end{align*}
The second integral in this expression is the normal derivative of a
constant function, and therefore is zero. The remaining integral has an
integrand with a singularity that also behaves as $\|{\mathbf{x}} - {\mathbf{y}}\|^{-1}$,
and odd-even integration can be applied.

Having developed a quadrature method to compute $\widehat{J}^H({\mathbf{x}};s)$,
the point-wise flux can be computed at time $t$ by applying the midpoint
rule along the Talbot contour in Fig.~\ref{eq:Talbot}. Then, the total
flux into ${\partial}\Omega$ can easily be computed by applying the trapezoid
rule to
\begin{align*}
  S(t) = \int_{{\partial}\Omega} {\mathcal{J}}(s,t) \, ds.
\end{align*}

\section{Particle based Kinetic Monte Carlo simulations}\label{sec:KMC}

Monte Carlo simulations provide a valuable tool for numerically estimating the distribution of capture times of diffusing particles for problems such as \eqref{eqn:IntroP}  and have been used extensively \cite{berez2014, Northrup1988,Litwin1980,Opplestrup2006}. Trajectories associated to the density \eqref{eqn:IntroP} can be constructed through the discretization
\begin{equation}\label{iter1}
\mathbf{x}(t+ dt) = \mathbf{x}(t) + \sqrt{2D dt} \, Z, \qquad \mathbf{x}(0) = \mathbf{x}_0,
\end{equation}
where $Z\sim N(0,1)$. The sequence of small displacements \eqref{iter1} terminates when the particle encounters the absorbing surface $\partial\Omega_a$. The algorithm is repeated for many particles (millions or even billions) to sample the capture time distribution. This approach is flexible and easy to implement but hampered by a set of problems.

If a fixed stepsize $dt$ is adopted, errors are introduced at that lengthscale that accrue near boundaries. First, for a capture event in the interval $(t,t+dt)$, we typically choose $t+dt$ as the arrival time which is an overestimate. Second, trajectories drawn from \eqref{iter1} will necessarily miss some encounters with boundaries and therefore overestimate the hitting time. Another challenge is that capture problems associated with \eqref{eqn:IntroP} are notorious for their fat-tailed distributions, i.e., a significant fraction of realizations undergo long excursions before capture. A key component of any
efficient method is adaptivity in stepsize since a trajectory of \eqref{iter1} simulated with a fixed step method will take a very long time to reach an absorbing site.

\subsection{Kinetic Monte Carlo (KMC) method for simulation of planar diffusion to absorbers}

Decreasing the step size in an adaptive manner based on distance to
target can ameliorate these issues. The Kinetic Monte Carlo (KMC) method
\cite{BatBerShv2003} maximizes this opportunity by advancing the
diffusion process in a spatial stepsize corresponding to the distance to
the target, $d(\mathbf{x}_0,\partial\Omega)$. The geometry of each step can take
many forms, but it should be chosen such that the details of the sojourn
can be rapidly and accurately sampled from closed form expressions.
Similar ideas have been employed in $N$-body simulations of kinetic
gases \cite{Opplestrup2006} and chemical reactions \cite{Wu2006}. In
this paper we describe implementation details and rudimentary analysis
of such a scheme that can handle complex geometries and mixed boundary
conditions. This method completely bypasses the need to advance
particles based on discretized steps such as~\eqref{iter1}.

\paragraph{Setup} We adopt a piecewise linear representation for the
boundary $\partial\Omega$ of the target set $\Omega$ based on vertexes
with $M$ straight edges $\partial\Omega_j$ so that $\partial \Omega =
\bigcup_{j=1}^M\partial\Omega_j$. On each boundary edge
$\partial\Omega_j$, we precalculate midpoints, unit normal vectors and
associate either Neumann or Dirichlet boundary conditions (others such
as Robin can be incorporated too). In addition, we calculate $R_0$, the
radius of the smallest circle centered at the origin that encloses all
targets (see Fig.~\ref{fig:WoS1b}).

A frequent and potentially expensive operation is the determination of
$d(\mathbf{x}_0,\partial\Omega)$, the distance of
$\mathbf{x}_0\in\RR^2\setminus\Omega$ to the nearest target. A simple approach
is to calculate the distance of $\mathbf{x}_0$ to each vertex of
$\partial\Omega$ and select the minimum. However, for highly refined
target geometries or numerous targets, the number of vertexes to scan
over may be prohibitive. To accommodate such scenarios, we employ a
particle-in-cell method that consists of a hierarchy of Cartesian grids
that envelop $\partial\Omega$ (Fig.~\ref{fig:quadtree}). The method
first queries the midpoints of the coarsest grid and uses simple
geometric criteria to eliminate those that cannot contain the closest
point. Queries are made of remaining subgrids at the next level of
refinement until a predefined level of refinement or a single point
remains. This results in a vastly smaller set of candidate vertexes to
calculate pointwise distances at the cost of some overheard and extends
this approach to large and complex targets sets.

\begin{figure}[htbp]
\centering
\includegraphics[width = 0.5\textwidth]{QuadTreeSchematic.pdf}
\caption{Visualization of the tree based particle-in-cell method for
  efficient evaluation of the distance $d(\mathbf{x}_0,\partial\Omega)$. The
  coloring of each point indicates the highest grid level that the
  vertex of $\partial\Omega$ is considered for the closest point
  evaluation. If the maximum refinement level is $k$, then all points in
  levels $1,2,\ldots, k-1$ are excluded from the distance calculation of
  $d(\mathbf{x}_0,\partial\Omega)$.
  \label{fig:quadtree}}
\end{figure}

With this setup in place, for each free particle $\mathbf{x}_0\in
\RR^2\setminus\Omega$, we calculate the shortest distance $R =
d(\mathbf{x}_0,\partial\Omega)$ and the associated projection $\textbf{p} =
\textrm{proj}_{\mathcal{P}_j} \mathbf{x}_0$ where $\mathcal{P}_j$ is the line that
contains the closest edge $\partial\Omega_j$. The position of the
particle is advanced based on four projection steps
(Fig.~\ref{fig:WoS1}) described below.
\begin{figure}[htbp]
\centering
\subfigure[Stage I: walk on sphere step.]{\includegraphics[width=0.3\textwidth]{2DWOSFIG1.pdf}\label{fig:WoS1a}}\qquad
\subfigure[Stage II: reinsertion step.]{\includegraphics[width=0.4\textwidth]{2DWOSFIG2.pdf}\label{fig:WoS1b}}
\subfigure[Stage III: square projection step.]{\includegraphics[width=0.4\textwidth]{SquareProjector.pdf}\label{fig:WoS1c}}\qquad
\subfigure[Stage IV: reflection step.]{\includegraphics[width=0.4\textwidth]{SemiCircleProjector.pdf}\label{fig:WoS1d}}
\caption{Schematic of the four stages of the planar KMC algorithm. \label{fig:WoS1}}
\end{figure}

\paragraph{Stage I: Radially symmetric projector}  If $R\in (R_{\textrm{min}},R_{\textrm{max}})$ such that the particle is neither too close nor too far from a target, we project to a ball of radius $R$ centered at $\mathbf{x}_0$. The parameters $R_{\textrm{min}},R_{\textrm{max}}$ are associated with stages II and III respectively and defined shortly. The time duration of this projection step is determined from the solution of the radial diffusion equation with a zero Dirichlet boundary condition at $r=R$ and a Dirac initial condition specifying the particle is initially at the origin. The solution $u(r,t)$ of the parabolic equation
\begin{align*}
\frac{\partial u}{\partial t} = D \frac{1}{r} \frac{\partial}{\partial
  r}\left(r \frac{\partial u}{\partial r} \right), \quad r\in(0,R), \quad t>0; \qquad u(R,t) = 0, \quad t>0; \qquad u(r,0) = \frac{\delta (r)}{r}, \quad r\in(0,R),
\end{align*}
gives the cumulative distribution of arrival times at $r=R$ to be
\begin{equation}\label{eq:2DArrivalTimes}
F(\tau) = 1- 2\sum_{n=0}^{\infty}\frac{ e^{-z_n^2 \tau}}{z_n J_1(z_n)}, \qquad t = \frac{D}{R^2} \tau, \qquad J_0(z_n) = 0, \quad n =0,1,2,\ldots.
\end{equation}
This distribution is sampled by drawing a uniform number $U\in(0,1)$ and solving $F(\tau) = U$. The CDF $F(\tau)$ is efficiently sampled by precomputing values of $z_n$ and $J_1(z_n)$ and using only as many terms as is necessary to approximate $F(\tau)$ to a predetermined tolerance.

\paragraph{Stage I (alternative for convex shapes): Plane projector} This projector can be adopted when the target is strictly convex so that the entire absorbing target lies to one side of the tangent plane to the surface. After translating the projection point ${\textbf{p}}$ to the origin and rotating by the slope of the incident edge, the projection step arises from the solution to the heat equation in the upper half plane with initial location $(0,y_0)$. Combining the relevant fundamental solutions with method of images yields the density
\begin{align*}
u(x,y,t) = \frac{1}{4\pi D t} \left[e^{- \frac{x^2 + (y-y_0)^2}{4Dt}} - e^{- \frac{x^2 + (y+y_0)^2}{4Dt}} \right],
\end{align*}
with association arrival time distribution to the plane
\begin{align*}
  \rho_T(t) = \int_{\RR} u_y(x,0,t) \,\mathrm{d}x =
    \frac{y_0}{2\sqrt{\pi} (Dt)^{\frac32}} e^{-\frac{y_0^2}{4Dt}}.
\end{align*}
The cumulative distribution is $\int_0^t \rho_T(\tau) \,\mathrm{d}\tau =
\text{erfc}(y_0/\sqrt{4Dt})$ so that the arrival time is sampled as
\begin{subequations}\label{eq:Plane}
\begin{align}\label{eq:PlaneTime}
t_{\ast} = \frac{1}{4D} \left[ \frac{y_0}{ \text{erfc}^{-1}(\eta)} \right]^2, \qquad \eta\in(0,1).
\end{align}
The hitting location on the tangent line is determined by the displacement $x_{\ast}$ from the projection point ${\bf p}$ which has the (Gaussian) distribution
\begin{align}\label{eq:PlaneLocation}
P_{X}(x_{\ast}) = \frac{u_y(x_{\ast},0,t_{\ast}) }{\rho_T(t_{\ast})} = \frac{1}{\sqrt{4\pi D t_{\ast}}} e^{-\frac{x_{\ast}^2}{4Dt_{\ast}} },
\end{align}
\end{subequations}$\!$
so that $x_{\ast} \sim \mathcal{N}(0,2Dt_{\ast})$.

\paragraph{Stage II: Reinsertion} It is inefficient to simulate the detailed trajectory of particles far from the absorbers. Therefore, when the distance $R = d(\mathbf{x}_0,\partial\Omega)$ exceeds a threshold ($R>  R_{\textrm{max}}$), we project the particle to a smaller disc of radius $R_{\textrm{ins}} > R_0$ that encloses all the absorbers (Fig.~\ref{fig:WoS1b}). Similar reinsertion procedures have been utilized in Monte Carlo solutions of elliptic problems \cite{OKMascagni2004,DOBRAMYSL2018}. Here we must take additional care to sample both the reinsertion time and the time dependent reinsertion location correctly. The arrival distribution for a particle initially on the $x$-axis (see Appendix \ref{sec:PlanarDist}) is
\begin{subequations}\label{flux:MAINomega}
\begin{align}\label{flux:MAINomega_a}
\mathcal{J}(\tau,\theta) =   \frac{D}{R_{\textrm{ins}}^2} \left[\frac{1}{2\pi} \chi_0(\tau)  + \frac{1}{\pi}\sum_{n=1}^{\infty}\chi_n(\tau) \cos n\theta\right] , \qquad \theta \in (-\pi,\pi), \quad  t>0, \quad \tau = \frac{D}{R_{\textrm{ins}}^2} t,
\end{align}
where the coefficients are
\begin{align}
  \label{flux:MAINomega_b}
  \chi_n(\tau) = \frac{2}{\pi}\int_{0}^{\infty}\left[
    \frac{J_n(\omega)Y_n(\omega \mu)-J_n(\omega\mu)Y_n(\omega)}
    {Y_{n}^{2}(\omega)+J_{n}^{2}(\omega)}\right]\omega
    e^{-\omega^2 \tau}\,d\omega.
\end{align}
\end{subequations}$\!$
Here $\mu := |\mathbf{x}_0|/R_{\textrm{ins}}$ is the ratio of the distance
$|\mathbf{x}_0|$ of the initial location from the origin to the reinsertion
radius $R_{\textrm{ins}}$. The optimal reinsertion radius is
$R_{\textrm{ins}} = R_0$ where $R_0$ is the radius of the smallest disc
enclosing all absorbers as shown in Fig.~\ref{fig:WoS1b}. However, many
computational efficiencies are gained by sampling \eqref{flux:MAINomega}
for a fixed value of $\mu$ - in practice we take $\mu = 60$,
$R_{\textrm{max}} = \mu R_0$ and reinsert to the disc of radius
$R_{\textrm{ins}} = |\mathbf{x}_0|/\mu$. By fixing $\mu$, the integrands of
\eqref{flux:MAINomega_b} can be tabulated over a range of $\omega$
values for efficient quadrature.

The first step in the sampling of \eqref{flux:MAINomega} is to determine
the arrival time density $\int_{\theta=0}^{2\pi} \mathcal{J}(
\tau,\theta)\,d\theta = \frac{D}{R_{\textrm{ins}}^2}\chi_0( \tau) $ with
associated CDF $F_T(t) = \int_0^{ \tau} \chi_0(\eta)\,d\eta$ where $\tau
= \frac{D}{R_{\textrm{ins}}^2}t$. The arrival time is sampled first
since the location will be dependent on this value---for shorter times
the arrival location is more tightly focussed around the initial
location while for larger arrival times, the insertion location has a
weaker dependence on the start location and approaches a uniform
distribution (see Fig.~\ref{fig:2DarrivalCDF_b}). For smaller values of
$\tau$\, (in practice $\tau<10^{10}$), the values of the integrand \eqref{flux:MAINomega_b} and
the associated CDF are tabulated over a range of $\omega$ values for
rapid quadrature. The sampling of $\chi_0(\tau)$ can be quite delicate
for large $\tau$ due to the slow convergence of the integral. To see
this, consider that for $\tau\gg1$ the main contribution to the integral
$\chi_0(\tau)$ is when $\omega^2 \tau = \mathcal{O}(1)$ or $\omega\ll1$. In
this regime we have that,
\begin{align}\label{eq:ApproxLarget}
\chi_0(\tau) \sim \frac{4 \log \mu}{\pi^2} \int_{0}^{\infty} \frac{w}{1+
  \frac{4}{\pi^2}(\gamma_e + \log(w/2))^2} e^{-\omega^2 \tau} \, d\omega
  = \mathcal{O}\left(\frac{1}{\tau|\log \tau|^2}\right), \qquad \tau \to \infty.
\end{align}
When $\tau\gg1$ (in practice $\tau >10^{10}$), we use the limiting form \eqref{eq:ApproxLarget} to posit the following explicit form of the density
\begin{align*}
  \chi_0(\tau) = \frac{4\log \mu}{\pi^2}
  \frac{1}{\tau(a_1 + a_2 \log \tau + a_3 \log \tau ^2)} +
  \mathcal{O}\left(\frac{1}{\tau^2}\right), \qquad \tau\gg1,
\end{align*}
for constants $a_1, a_2, a_3$ determined from fitting. This gives the exact cumulative density function for $\tau\gg1$
\begin{align}\label{eq:fittedDensity}
F_T(t) = 1- \int_{\tau}^{\infty} \chi_0(\eta)\,d\eta = 1 - \frac{4\log
  \mu}{\pi^2\sqrt{4a_1 a_3 - a_2^2}} \left( \pi - 2\tan^{-1}
  \left[\frac{a_2 + 2a_3 \log \tau}{\sqrt{4a_1 a_3 - a_2^2}}\right]
  \right), \quad \tau = \frac{D}{R_{\textrm{ins}}^2} t,
\end{align}
where for $\mu = 60$, we obtain from fitting the constants
\begin{align*}
a_1 = 1.4670, \qquad a_2 = 0.3102, \qquad a_3 = 0.2029.
\end{align*}
\begin{figure}[htbp]
\centering
\subfigure[Arrival time distribution.]{\includegraphics[width = 0.4\textwidth]{2DArrivalTimeCDF.pdf}\label{fig:2DarrivalCDF_a}}\hspace{0.5in}
\subfigure[Arrival angle distribution for several $t_{\ast}$ values.]{\includegraphics[width = 0.39\textwidth]{2DArrivalAngleCDF.pdf}\label{fig:2DarrivalCDF_b}}
\caption{CDFs of arrival time and location for the reinsertion step
  when $D=1$, $\mu=60$. Left: The reinsertion time distribution
  $F_{\theta}(\tau) = \int_0^{\tau} \chi_0(\eta)\,d\eta$ given by
  \eqref{eq:dist2Dreinsertion} (dashed line indicates the cutoff
  $\tau =10^{10}$ after which the fitted density
  \eqref{eq:fittedDensity} is used). Right: The arrival location CDF
  \eqref{eq:samp:Th} for several arrival times $\tau_{\ast}$. The distribution is largely uniform for large $\tau_{\ast}$. \label{fig:2DarrivalCDF}}
\end{figure}
For a particular arrival time realization $\tau = \tau_{\ast}$, the angular location $\theta$ of reinsertion satisfies
\begin{equation}\label{eq:samp:Th}
F_{\theta}(\theta; \tau_{\ast}) = \int_{-\pi}^{\theta}
  \frac{\mathcal{J}(\tau_{\ast},\eta)}{\chi_0(\tau_{\ast})}\,d\eta = \frac{\theta+\pi}{2\pi} + \frac{1}{\pi}\sum_{n=1}^{\infty} \frac{\sin (n \theta)}{n} \frac{\chi_n( \tau_{\ast})}{\chi_0(\tau_{\ast})}, \qquad \theta\in(-\pi,\pi).
\end{equation}
To sample from \eqref{eq:samp:Th}, a uniform number $U\in(0,1)$ is drawn and the equation $F_{\theta}(\theta;\tau_{\ast}) = U$ solved for $\theta$. In practice we use an adaptive procedure where we retain only the terms satisfying $|\chi_n(\tau_{\ast})/(n\chi_0(\tau_{\ast}))| > \tau_{\textrm{tol}}$ in the summation of \eqref{eq:samp:Th}. For a large proportion of realizations, the arrival time $\tau_{\ast}$ is sufficiently large (in practice $\tau_{\ast}>10^{10}$) so that the first term is negligible and the arrival location is uniformly distributed on the disc. In Fig.~\ref{fig:2DarrivalCDF}, we plot the CDFs of arrival time distribution and arrival location distribution for $\mu=60$. This method permits rapid and accurate sampling of the reinsertion step.

\paragraph{Stage III: Square projector}
If $R<R_{\textrm{min}}$, then the particle is close enough to determine if contact occurs. By \lq\lq close enough\rq\rq, we mean that the projection $\bf{p}$ lies within the edge segment so that a square of side length $2R$ centered at $\mathbf{x}_0$ lies entirely within the target edge  (cf.~Fig.~\ref{fig:WoS1c}). This gives explicitly that $R_{\textrm{min}} = \min(d_1,d_2)$ where $d_1$, $d_2$ are the distances between $\bf{p}$ and the edge vertexes (cf.~Fig.~\ref{fig:WoS1c}). The projection step is then determined from the solution to the parabolic equation on the square $\mathcal{S} = [-R, R]^2$
\begin{subequations}\label{eq:SquareDiffusion2D}
\begin{gather}
 \label{eq:SquareDiffusion2D_a} \frac{\partial u}{\partial t} = D \left(
  \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial
  y^2}\right), \quad \mathbf{x} \in \mathcal{S}, \quad t > 0; \quad u(\mathbf{x}, t) = 0, \quad \mathbf{x} \in \partial \mathcal{S}, \quad t>0; \\[5pt] \quad u(\mathbf{x}, 0) = \delta(x)\delta(y), \quad \mathbf{x} \in \mathcal{S}, \quad t= 0.
\end{gather}
\end{subequations}$\!$
The separable solution to \eqref{eq:SquareDiffusion2D} yields the CDF of first arrival times to the square edge $\partial\mathcal{S}$
\begin{align*}
P_T(\tau) = \int_0^{\tau} \rho_T(\eta)\,d\eta = 1-\frac{32}{\pi^2}  \sum_{l=0}^{\infty} \sum_{k=0}^{\infty}  \frac{2k+1}{2l+1} \frac{(-1)^{l+k}}{(2k+1)^2 + (2l+1)^2} e^{-((2l+1)^2 + (2k+1)^2) \pi^2 \tau}. \quad \tau = \frac{D}{R^2}t.
\end{align*}
This distribution is sampled by drawing a uniform number $U\in(0,1)$ and solving $P_T(\tau^{\ast}) = U$. Each side of the square has an equal probability $1/4$ of being hit and the arrival location is sampled from the density
  \begin{align*}
\rho_X(x) = \frac{4u_y(x,0,\tau^{\ast})}{\rho_T(\tau^{\ast})} = \frac{\pi}{2}\, \ds\frac{\ds\sum_{l=0}^{\infty} \ds\sum_{k=0}^{\infty}  (2l+1) (-1)^{l+k}  \, e^{-((2l+1)^2 + (2k+1)^2) \pi^2 \tau^{\ast}}\sin[(2k+1)\pi x]}{\ds\sum_{l=0}^{\infty} \ds\sum_{k=0}^{\infty}  \frac{2k+1}{2l+1} (-1)^{l+k} e^{-((2l+1)^2 + (2k+1)^2) \pi^2 \tau^{\ast}}}.
\end{align*}
In practice, we precalculate a large number of pairs $\{(\tau_j,x_j)\}$ to draw from when needed during simulation.
\paragraph{Stage IV: Reflection step}

In the scenario that the particle hits a reflecting portion of the target surface, it is then projected back into the bulk onto a semi-circle of radius $R = \min(d_1,d_2)$, corresponding to the distance to the nearest vertex. In practice, we avoid rounding errors by setting $R = \max(\epsilon,\min(d_1,d_2))$ where $\epsilon$ is a small number comparable to machine precision. In the reflecting boundary condition scenario, the projection step is identical to that of stage I with uniform location and arrival time sampled from \eqref{eq:2DArrivalTimes}. For highly convoluted geometries, it may be that this semi-circle intersects with distal elements of the target. This scenario can be accounted for by calculating the shortest distance $d$ to other segments of the boundary and setting $R = \max(\epsilon,\min(d_1,d_2,d))$.

\subsection{Rudimentary convergence analysis in half plane case}

In this section we present some analysis of the convergence properties of the KMC approach in the simplified scenario of diffusion in the upper half plane with capture in the window $\{(x,0) \, \mid \,  |x|<h/2\}$ with particular emphasis on the role of reinsertion. Specifically, we solve the equation
\begin{subequations}\label{eq:half_plane_full}
\begin{gather}
\label{eq:half_plane_full_a}    \frac{\partial u}{\partial t} =
  \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2},
  \quad x\in\RR, \quad y >0, \quad t>0; \\[5pt]
\label{eq:half_plane_full_b}     u = 1, \quad y = 0, \quad|x| < \frac{h}{2}, \qquad u_y = 0,\quad y = 0, \quad |x| \geq \frac{h}{2};\\[5pt]
\label{eq:half_plane_full_c}  u = \delta(x-x_0)\delta(y-y_0), \quad
  x\in\RR, \quad y >0, \quad t=0.
\end{gather}
\end{subequations}$\!$
with a simplified KMC method composed of Stages 1(alternative) and 3.
From an initial point $\mathbf{x}_0=(x_0,y_0)$, this algorithm results in a
sequence of points $(\mathbf{x}_1,\mathbf{x}_2,\ldots, \mathbf{x}_n)\in \RR^{2+}$ which will eventually alight on the absorbing portion.  For an ensemble of particles, we denote $a_n$ to be the fraction free after $n$ iterations so
\[
a_{n+1} =  (1-{p_n}) a_{n}, \qquad a_0 = 1,
\]
where ${p_n}$ is the probability of capture at the $n^{th}$ iteration. To investigate $p_n$, we first consider the splitting problem $p(x,y)$ for the probability that a particle starting at $(x,y)$ first contacts the plane $y=0$ on the absorbing window. This satisfies
\begin{subequations}\label{split}
\begin{gather}
\label{split_a}    \frac{\partial^2 p}{\partial x^2} +  \frac{\partial^2
  p}{\partial y^2} = 0, \quad x\in\RR, \quad y>0;\\[5pt]
\label{split_b}       p = 1, \qquad y = 0, \quad |x|< h/2, \qquad p =0, \qquad y =0, \quad |x|\geq h/2,
\end{gather}
and admits the solution
\begin{equation}\label{split_c}
    p(x,y) = \frac{1}{\pi}\left[ \tan^{-1}\left(\frac{x+\frac{h}{2}}{y} \right) - \tan^{-1}\left(\frac{x-\frac{h}{2}}{y} \right)  \right].
\end{equation}
\end{subequations}$\!$
We consider that the current location $\mathbf{x}_n = (x_n,y_n)$ has arisen from a projection step (stage III). Without loss of generality, we assume the previous contact $(\bar{x}_n,0)$ with the $y_n=0$ plane is such that $\bar{x}_n> \frac{h}{2}$. We may then parameterize (cf.~Fig.~\ref{Fig:reinsertion_schematic}) the point $\mathbf{x}_n = (x_n,y_n)$ as
\begin{align*}
x_n = \frac{h}{2} + \bar{r}_n(1+\cos\phi), \qquad y_n = \bar{r}_n\sin\phi, \qquad \phi \in(0,\pi), \qquad \bar{x}\in(h/2,\infty).
\end{align*}

\begin{figure}[H]
\centering
\subfigure[KMC step]{\includegraphics[clip, width = 0.425\textwidth]{ReinsertionSchematic.pdf}}
\qquad
\subfigure[Reinsertion]{\includegraphics[clip, width = 0.45\textwidth]{HalfPlaneReinsertion.png}}
\caption{A schematic showing the simplified KMC solution in the upper half place. (a) the projection step from the real line back into the bulk. (b) The reinsertion step to keep the particle in the vicinity of the absorbing set. \label{Fig:reinsertion_schematic}}
\end{figure}
It follows from the splitting probability \eqref{split_c} that
\begin{align}
\nonumber p(x_n,y_n) &= \frac{1}{\pi}\left[ \tan^{-1} \left( \frac{x_n-\frac{h}{2}}{y_n} + \frac{h}{y_n} \right) - \tan^{-1} \left( \frac{x_n-\frac{h}{2}}{y_n}\right)\right]\\[5pt]
\nonumber {} & = \frac{1}{\pi} \left[  \tan^{-1} \left( \cot \frac{\phi}{2} + \frac{h}{\bar{r}_n} \csc \phi \right) - \tan^{-1} \left( \cot \frac{\phi}{2}\right) \right]\\[5pt]
\nonumber {} & = \frac{1}{\pi} \left[  \tan^{-1} \left( \tan \left(\frac{\pi}{2} - \frac{\phi}{2}\right) + \frac{h}{\bar{r}_n} \csc \phi \right) - \left(\frac{\pi}{2} - \frac{\phi}{2}\right) \right],
\end{align}
where the angles $\phi\in(0,\pi)$ are distributed uniformly. Applying the change of variables $\eta = \frac{\pi-\phi}{2}$ yields
\begin{equation*}
    p = \frac{1}{\pi} \left[  \tan^{-1} \left( \tan \eta + \alpha_n \csc 2\eta \right) - \eta \right], \qquad \alpha_n = \frac{h}{\bar{r}_n}.
\end{equation*}
We now define the average probability $p_n =
\frac{1}{\pi}\int_{\phi=0}^{\pi} p \, d\phi =
\frac{2}{\pi}\int_{\eta=0}^{\pi/2} p \, d\eta$ and find that
\begin{align}\label{eqn:pn}
p_n = \frac{2}{\pi^2}\int_{\eta = 0}^{\frac{\pi}{2}} \left[ \tan^{-1}
  \left( \tan \eta + \alpha_n \csc 2\eta \right) - \eta \right] \, d\eta, \qquad \alpha_n = \frac{h}{\bar{r}_n}.
\end{align}

Without reinsertion, many trajectories will yield very small values for the parameter $\alpha_n$ as $\bar{r}_n$ can attain very large values. To gain further insight into this scenario, we consider the limiting case of \eqref{eqn:pn} as $\alpha\to0$.

\subsection{Asymptotic analysis of splitting probability}

Here we develop an asymptotic approximation for \eqref{eqn:pn} in the limit as $\alpha\to0$ (subscript $n$ dropped for convenience). The integral features global contributions and local contributions near $\eta=0$. To delineate between these contributions, we define the small parameter $\delta>0$ such that $\alpha \ll \delta \ll 1$. Then we have that
\begin{align}\label{eqn:pn_alpha}
    p_n(\alpha) = \underbrace{\frac{2}{\pi^2}\int_{\eta = 0}^{\delta}
    \left[ \tan^{-1} \left( \tan \eta + \alpha \csc 2\eta \right) - \eta
    \right] \, d\eta}_{A_1} +
    \underbrace{\frac{2}{\pi^2}\int_{\eta = \delta}^{\frac{\pi}{2}}
    \left[ \tan^{-1} \left( \tan \eta + \alpha \csc 2\eta \right) - \eta
    \right] \, d\eta}_{A_2}.
\end{align}
where $A_1$ and $A_2$ will be considered separately and combined so that their sum is independent of $\delta$.

\paragraph{Evaluation of $A_1$} In this region we apply the
approximation $\tan^{-1}(x + y) = \tan^{-1}(x) + y/(1+x^2) + \mathcal{O}(y^2)$ for $y\to0$ with $x = \alpha \csc 2\eta$ and $y = \tan \eta$. Then we have that
\begin{equation}
 \nonumber A_1 := \frac{2}{\pi^2}\int_{0}^{\delta} \left[ \tan^{-1}
  \left( \tan \eta + \alpha \csc 2\eta \right) - \eta \right]\, d\eta \sim
  \frac{2}{\pi^2} \int_{0}^{\delta} \left[ \tan^{-1} (\alpha\csc 2\eta)
  + \frac{\tan \eta}{1+ \alpha^2 \csc^2(2\eta)} - \eta\right]\, d\eta.
\end{equation}
Applying small argument approximations for $\eta\ll1$, we have that
\begin{align}
\nonumber A_1 &\sim \frac{2}{\pi^2}\int_{0}^{\delta}\left(
  \tan^{-1}\left(\frac{\alpha}{2\eta}\right)  - \frac{\alpha^2
  \eta}{4\eta^2 + \alpha^2}\right) \, d \eta = \frac{\alpha}{\pi^2}
  \int_0^{\tfrac{2\delta}{\alpha}}\left( \tan^{-1} \frac{1}{z} -
  \frac{\alpha}{2} \frac{z}{z^2+1}\right) \, dz, \qquad (\eta = \alpha z/2)\\[5pt]
\nonumber {} & =  \frac{\alpha}{\pi^2} \left[\frac{1}{2} \log(1+ z^2) +
  z \tan^{-1}\frac{1}{z} - \frac{\alpha}{4}\log(1+ z^2)
  \right]_0^{\frac{2\delta}{\alpha}} = \frac{\alpha}{\pi^2} \left[
    \left(\frac{1}{2} - \frac{\alpha}{4}\right) \log\left(1 +
    \frac{4\delta^2}{\alpha^2}\right) + \frac{2\delta}{\alpha} \tan^{-1}
    \frac{\alpha}{2\delta}  \right].
\end{align}
The parameter $\delta$ is defined so that $\alpha \ll \delta$ and therefore $\frac{\delta}{\alpha}\gg1$. Applying Taylor we find that
\begin{align}\label{eqn:A1}
A_1 = \frac{\alpha}{\pi^2} \left( \log \frac{2\delta}{\alpha} + 1
  \right) + \mathcal{O}\left(\frac{\alpha^3}{\delta^2},
  \alpha^2\log\left(\frac{\delta}{{\displaystyle \varepsilon}}\right)\right).
\end{align}

\paragraph{Evaluation of $A_2$} In the integral $A_2$ we
have that $\alpha\ll\delta$ and apply the approximation $\tan^{-1}(x +
y) = \tan^{-1}(x) + y/(1+x^2) + \mathcal{O}(y^2)$ for $y\to0$ with $x = \tan \eta$ and $y = \alpha \csc 2\eta$. It then follows that
\begin{align}
\nonumber A_2 &:= \frac{2}{\pi^2}\int_{\eta = \delta}^{\frac{\pi}{2}}
  \left[ \tan^{-1} \left( \tan \eta + \alpha \csc 2\eta \right) - \eta
  \right] \, d\eta \sim \frac{2}{\pi^2}\int_{\delta}^{\frac{\pi}{2}}
  \left[\tan^{-1} (\tan\eta) + \alpha \frac{\csc 2\eta}{1 + \tan^2\eta}
  - \eta \right] \, d\eta\\[5pt]
\nonumber & = \frac{2\alpha}{\pi^2}\int_{\delta}^{\frac{\pi}{2}}
  \frac{\csc 2\eta}{1 + \tan^2\eta} \, d\eta = \frac{\alpha}{\pi^2}
  \int_{\delta}^{\frac{\pi}{2}} \cot\eta \, d\eta = \frac{\alpha}{\pi^2}
  \left[\log\sin\eta \right]_{\delta}^{\frac{\pi}{2}}\\[5pt]
& = -\frac{\alpha}{\pi^2} \log\delta. \label{eqn:A2}
\end{align}

To finalize the approximation of \eqref{eqn:pn_alpha}, we combine expressions \eqref{eqn:A1} and \eqref{eqn:A2} and reintroduce $\alpha_n$ reflecting that this parameter changes over each iteration. This yields that
\begin{equation}\label{eqn:pn_asy}
    p_n \sim A_1 + A_2 = \frac{\alpha_n}{\pi^2} \left(\log \frac{2}{\alpha_n} + 1 \right), \qquad \alpha_n \to 0, \qquad \alpha_n = \frac{h}{\bar{r}_n}.
\end{equation}
Hence we see that while the probability of capture is positive at each iteration, it can become arbitrarily small as $|\bar{x}_n| =  \frac{h}{2} + \bar{r}_n \to\infty$. This results in an algorithm with polynomial convergence rate (see Fig.~\ref{fig:reinsertion_convergence_a}).

\begin{figure}
    \centering
\subfigure[Capture probability without reinsertion.]{\includegraphics[width= 0.45\textwidth]{AsyCompProb.png}}
\subfigure[Capture probability from reinsertion disk.]{\includegraphics[width= 0.45\textwidth]{AsyCompProb2.png}}
    \caption{(Left figure) Comparison of the capture probability $p_n(\alpha)$ (solid curve) \eqref{eqn:pn} and the small arguement approximation \eqref{eqn:pn_asy} (dashed curve) without reinsertion. We note that $p\to 1/4$ for $\alpha_n = h/\bar{r}_n \to \infty$. (Right figure) Here we show the probability of capture with reinsertion to a disk of radius $R$. }
    \label{fig:my_label}
\end{figure}

\subsection{Reinsertion analysis}

The aim of reinsertion is to reestablish exponential convergence rate in the algorithm by limiting the maximum value of $\bar{x}_n$ and hence promoting faster capture.

Reinsertion projects wayward particles back to a smaller disc of radius $R$ that encloses all targets and therefore omits simulating  trajectories far from the capture regions. When reinserting from sufficiently distant points, the placement on the disc is largely uniform with $(x,y) = R(\cos\theta,\sin\theta)$ with small corrections given by \eqref{eq:samp:Th}. Points on this disk have average probability of capture (see \eqref{split})
\begin{align}
\nonumber p_n &=
   \frac{1}{\pi}\int_{0}^{\pi}p(R\cos\theta,R\sin\theta)\,d\theta =
  \frac{1}{\pi^2} \int_{0}^{\pi} \left[  \tan^{-1} \big( \cot \theta +
  \frac{h}{2R} \csc \theta \big) - \tan^{-1} \big( \cot \theta -
  \frac{h}{2R} \csc \theta \big) \right]\, d\theta\\[5pt]
\label{prob_Rins} &\sim \frac{2}{\pi^2} \frac{h}{R}, \quad \mbox{as} \quad \frac{h}{R} \to 0.
\end{align}
In an ensemble of particles, some will be reinserted to the disk $x^2+y^2 = R^2$ while others remain inside it. Those inside have greater probability of capture at the next step, therefore the quantity \eqref{prob_Rins} reflects a lower bound on the likelihood of capture. Hence we see that the probability of capture after $n$ stages has bound $p_n \gtrapprox \frac{2}{\pi^2} \frac{h}{R}$. The key observation here is that the probability of capture at each iteration is now bounded below by a constant defined in term of the geometric parameter $h$ and the reinsertion radius $R$.  This ensures an exponential convergence rate of the algorithm with slower rates associated with smaller targets ($h\ll1$) and a faster rate associated with a smaller reinsertion radius $R$.

As an exposition of this analysis, we show in Fig.~\ref{fig:reinsertion_convergence} the convergence of the KMC method on the simplified half plane problem \eqref{eq:half_plane_full} with and without reinsertion. In the absence of reinsertion, polynomial convergence is attained as shown by linear behavior on a log-log plot. When reinsertion is implemented (to radius $R = 1$), we observe exponential convergence shown by linear behavior on a log plot. This demonstrates the key role of reinsertion in attaining exponential convergence of the KMC method.
\begin{figure}[htbp]
\centering
\subfigure[Convergence without reinsertion]{\includegraphics[clip, width
  = 0.45\textwidth]{ConvergenceWithoutReinsertion.pdf}\label{fig:reinsertion_convergence_a}}
\qquad
\subfigure[Convergence with reinsertion]{\includegraphics[clip, width =
  0.45\textwidth]{ConvergencewReinsertion.pdf}\label{fig:reinsertion_convergence_b}}
\caption{Convergence statistics of the KMC method for the simplified problem \eqref{eq:half_plane_full}. Left: In the absence of reinsertion, a polynomial rate of convergence is attained shown by the linear behavior on the log-log plot. Right: Reinsertion recovers exponential convergence as seen in the linear behavior on the log plot. Simulations based on ensemble of $10^6$ particles. \label{fig:reinsertion_convergence}}
\end{figure}

\section{Results}\label{sec:results}

In this section we show a variety of examples to demonstrate the ability of these methodologies for approximating full arrival time distributions to complex absorbing sets. In our simulations, we have used a common diffusivity of $D=1$. Arrival times from the KMC method are translated with $t\to t+1$ so that $t=0$ is mapped to $10^0$ in log space. \texttt{MATLAB}'s histogram function is then applied with the \lq\lq probability\rq\rq\ normalization option. The hybrid approaches refer to solving the Laplace transform using either the expansion asymptotic \eqref{sysmain} or the boundary integral method (BIM) described in Sec.~\ref{sec:BIE}. This is followed by numerical inversion of the transform equation as described in Sec.~\ref{sec:invLT}.

\subsection{Planar results}\label{sec:resultsPlanar}

Here we use three examples to validate the numerical KMC method and corroborate with both the hybrid approaches (asymptotic and boundary integral). The first example is a simple one target scenario in which the closed form solution \eqref{eqn:timepdf} is available for comparison. The remaining two examples show the efficacy of the method on more complex absorbing sets consisting of multiple targets of varying radii. The parameter values for the three examples are
\begin{align*}
(\mbox{One target}) & \quad \mathbf{x}_1 = (0,0), \quad r_1 = 0.05, \quad \mathbf{x}_0 = (5,0).\\[5pt]
(\mbox{Three targets})& \quad \mathbf{x}_1 = (3,3), \quad \mathbf{x}_2 = (8,8), \quad \mathbf{x}_3 = (10,10), \qquad r_{1}= \frac18, \quad r_{2} = \frac16, \quad r_{3} = \frac13, \quad \mathbf{x}_0 = (0,0).\\[5pt]
(\mbox{Six targets})& \quad \mathbf{x}_1 = (-3,0), \quad  \mathbf{x}_2 =(0,-2), \quad  \mathbf{x}_3 =(\sqrt{2},-\sqrt{2}),\quad \mathbf{x}_4 =(2,0),\\[5pt]
\phantom{(\mbox{Six targets})}& \quad \mathbf{x}_5 =(\sqrt{2},\sqrt{2}),\quad \mathbf{x}_6 =(0,2), \qquad
r_1  = 0.275, \quad r_{2} = \ldots =r_{6} = 0.02, \quad \mathbf{x}_0 = (0,0).
\end{align*}

The results shown in Fig.~\ref{fig:PlanarExs} show good agreement between the three approaches. In Fig.~\ref{fig:PlanarExs}(a) we compare the hybrid-asymptotic, KMC and exact one-pore solutions showing excellent agreement. In the two more challenging examples, we generally see good agreement between the asymptotic and boundary integral approaches. In the more challenging 6 target case, we see in Fig.~\ref{fig:PlanarExs}(c) a slightly diminished agreement is observed near the peak.
\begin{figure}[htbp]
 \vspace{-1\baselineskip
  \centering
   \subfigure[One target example]{\includegraphics[width=0.325\linewidth]{OnePoreEx.pdf}}
   \subfigure[Three target example]{\includegraphics[width=0.325\linewidth]{3PoreEx.pdf}}
   \subfigure[Six target example]{\includegraphics[width=0.325\linewidth]{6PoreEx.pdf}}
      \caption{Relative frequency of arrival times to planar targets by KMC, exact solution (one target case), and hybrid approaches (asymptotic and boundary integral). Schematics of target arrangements shown as insets.\label{fig:PlanarExs} }
\end{figure}

\subsection{Homogenization}\label{sec:resultsHomogenization}

In this example, we consider the case of a single target with a mix of absorbing and reflecting portions. This describes the scenario where an impermeable cellular membrane surface is covered in surface receptors.  We determine the full distribution of arrival times using both the KMC method and the hybrid asymptotic-homogenization result \eqref{eqn:homogenization}. In the application of the hybrid approach, we use the analytically determined logarithmic capacitance \eqref{eqn:homogenization}, or equivalently the effective radius, in the single patch result.

\begin{figure}[htbp]
    \centering
    \subfigure[Ex 1: $M = 12$, $f = 1/3$.]{\includegraphics[width = 0.495\textwidth]{HomogInner1-3.png}}
    \subfigure[Ex 1: $M = 12$, $f = 1/3$.]{\includegraphics[width = 0.495\textwidth]{FullDist1-3.png}}\\
        \subfigure[Ex 2: $M = 8$, $f = 1/8$.]{\includegraphics[width = 0.495\textwidth]{HomogInner1-8.png}}
    \subfigure[Ex 2: $M = 8$, $f = 1/8$.]{\includegraphics[width = 0.495\textwidth]{FullDist1-8.png}}
    \caption{Results from the KMC and asymptotic method on the homogenization examples \eqref{ex:Homog}. Left panels: target (scaled to unity radius) with layout of absorbers. Contours indicate the numerical solution of \eqref{eqn:LogCap}. Right figures: Agreement between the full arrival time densities obtained from KMC and the hybrid-asymptotic method.}
    \label{fig:HomoG}
\end{figure}

In the two examples shown in Fig.~\ref{fig:HomoG}, we take a circular absorbing target centered at the origin with radius ${\displaystyle \varepsilon} = 0.05$. The target itself features $M$ equally spaced absorbing windows centered at the roots of unity $(\cos \frac{2\pi k}{M},\sin \frac{2\pi k}{M} )$. The windows occupy a combined fraction $f$ and each has common angular extent $\sigma = \frac{2\pi f}{M}$. In each case we use the homogenized formula \eqref{eqn:homogenization_c} to obtain the logarithmic capacitance $d$ and then apply the result \eqref{2DResults} for $N=1$.

The relevant parameters obtained for the two examples are
\begin{subequations}\label{ex:Homog}
\begin{gather}
\label{ex:Homog_a} \text{Ex 1:} \quad M = 12, \quad f = \frac{1}{3}, \quad d = 0.8978, \quad \mathbf{x}_0 = [5,0]; \\[5pt]
\label{ex:Homog_b} \text{Ex 2:} \quad M = 8,\quad f =  \frac{1}{8}, \quad  d = 0.6657, \quad \mathbf{x}_0 = [2,0].
\end{gather}
\end{subequations}$\!$
In Fig.~\ref{fig:HomoG}(b),(d) we show good agreement between the two methods while in Fig.~\ref{fig:HomoG}(a),(c) we display a visualization of the solution to the capacitance problem \eqref{eqn:LogCap}. We remark that the homogenization effect can be seen through the gradual radial symmetrization of the contours.

\subsection{Clustering}\label{sec:resultsClustering}

In this section, we use two examples of clustered target configurations to compare results from the KMC method with both the asymptotic and BIM approaches. In addition, we show the effectiveness of homogenization where the clustered target configuration is replaced by a single circular target of appropriately chosen radius. The specific parameters are given by
\begin{subequations}\label{ex:cluster}
\begin{align}
\text{Ex 1:} \quad \mathbf{x}_{k} &=  \frac{3}{2}\left(\cos \frac{2\pi k}{5},\sin \frac{2\pi
  k}{5}\right), \quad k = 1,\ldots,5, \quad \mathbf{x}_6 = (0,0), \qquad  r_{1,\ldots,6} = 0.1, \quad \mathbf{x}_0 = (10,0);  \\[5pt]
\text{Ex 2:} \quad \mathbf{x}_{k} &=  \frac{3}{2}\left(\cos \frac{2\pi k}{8},\sin \frac{2\pi
  k}{8}\right), \quad k = 1,\ldots,8, \qquad r_{1,\ldots,8} = 0.05, \quad \mathbf{x}_0 = (5,0).
\end{align}
\end{subequations}$\!$

\begin{figure}[htbp]
\centering
\subfigure[Numerical solution of \eqref{eqn:LogCap} yielding logarithmic capacitance $d = 1.2295$.]{\includegraphics[width = 0.23\textwidth]{5ClusterandCenterCap.pdf}}\quad
\subfigure[Normalized target fluxes.]{\includegraphics[width = 0.32\textwidth]{5ClusterandCenterfluxes.pdf}}\quad
\subfigure[Agreement of capture times (from KMC) with asymptotic and homogenized densities.]{\includegraphics[width = 0.35\textwidth]{5ClusterandCenter.pdf}}\\
\subfigure[Numerical solution of \eqref{eqn:LogCap} yielding logarithmic capacitance $d = 1.2737$.]{\includegraphics[width = 0.24\textwidth]{8ClusterandCenterCap.pdf}} \quad
\subfigure[Normalized target fluxes.]{\includegraphics[width = 0.32\textwidth]{8ClusterFluxes.pdf}}\quad
\subfigure[Capture times (from KMC) with asymptotic and homogenized densities.]{\includegraphics[width = 0.35\textwidth]{8Cluster.pdf}}
\caption{Comparison of the full asymptotic density with those arising from the KMC method, boundary integral method (BIM) and homogenization (replacing the target cluster with a single absorbing target). }
\label{fig:Clustering}
\end{figure}
The numerical method described by equation \eqref{eqn:complexSeries} allows for the computation of the logarithmic capacitance for each of the absorbing sets (cf.~Fig.~\ref{fig:Clustering}(a,d)). For the parameters specified in \eqref{ex:cluster}, we determine $d = 1.2295$ (\text{Ex 1}) and $d= 1.2737$ (\text{Ex 2}). The logarithmic capacitance can be interpreted as the effective radius of a single target that reflects the capture potential of the cluster. As with the homogenization example in Sec \ref{sec:resultsHomogenization}, we observe that replacing complex configurations with a single target of appropriately chosen radius produces a very accurate representation of the full arrival time distribution. However, the single target representation does reduce certain direction information encoded in the distribution of arrivals over the targets in the cluster. The hybrid-asymptotic method can rapidly determine the fluxes by numerical inverse Laplace transform of $\hat{\mathcal{J}}_j =  2\pi D \nu_j S_j$  where each $S_j$ is determined by \eqref{sysmain}. In Fig.~\ref{fig:Clustering}(b,e) we show normalized fluxes into four sets of traps that are symmetrically arranged with respect to the initial location. The traps aligned towards the initial data accrue most of the inbound flux while the peaks, representing most likely arrival times at a particular target, are ordered by their distance to the initial location. The distribution of fluxes over the targets encode directional information that can infer the source location \cite{Miles2020}.

\subsection{Splitting Probabilities}\label{sec:resultsSplitting}

In this section, we demonstrate the convergence of the dynamic fluxes to the static splitting probabilities $\{\phi_k(\mathbf{x})\}_{k=1}^N$, where $\phi_k(\mathbf{x})$ denotes the probability that a diffusing particle originally at $\mathbf{x}$ reaches the $k^{th}$ target before any others. For exposition purposes, we focus on the scenario of completely absorbing targets. These probabilities satisfy the exterior Laplace problem
\begin{subequations}\label{eq:splits}
\begin{gather}
\label{eq:splits_a} \Delta \phi_k = 0, \quad \mathbf{x}\in\mathbb{R}^2\setminus\Omega; \qquad \phi_k(\mathbf{x}) = \delta_{jk},\quad \mathbf{x} \in\partial\mathcal{A}_j, \qquad j= 1,\ldots, N;\\[5pt]
\label{eq:splits_b} \phi_{k}(\mathbf{x}) \quad \text{finite as} \quad |\mathbf{x}| \to \infty,
\end{gather}
\end{subequations}$\!$
where $\delta_{jk}$ is the Kronecker delta. The asymptotic solution of \eqref{eq:splits} as ${\displaystyle \varepsilon}\to0$ is developed along similar lines to Sec.~\ref{sec:asy} (also see \cite[Sec 5]{Kurella2015}). Accordingly, we present the solution as ${\displaystyle \varepsilon}\to0$ directly as
\begin{subequations}\label{eqn:splitsol}
\begin{equation}\label{eqn:splitsol1}
\phi_k(\mathbf{x}) \sim 2\pi \sum_{j=1}^{N} \nu_j A_j G_0(\mathbf{x};\mathbf{x}_j) + \bar{\phi}, \qquad G_0(\mathbf{x};\mathbf{y}) = \frac{1}{2\pi}\log|\mathbf{x}-\mathbf{y}|,
\end{equation}\label{eqn:splitsol2}
where the $(N+1)$ constants $(A_1,\ldots,A_N, \bar{\phi})$ are determined from the linear system,
\begin{equation}
\sum_{j=1}^{N} \nu_j A_j = 0; \qquad -A_j + 2\pi \sum_{\substack{i=1\\ i\neq j}}^N \nu_i A_i G_0(\mathbf{x}_j;\mathbf{x}_i)  + \bar{\phi} = \delta_{jk}, \qquad j = 1,\ldots,N.
\end{equation}
\end{subequations}$\!$
The gauge functions $\nu_{j} = -1/\log{\displaystyle \varepsilon} d_j$ are defined in \eqref{eq:guage}. We remark that since capture is guaranteed for planar Brownian motion, we have that $\sum_{k=1}^N \phi_k(\mathbf{x}) = 1$.

\begin{figure}[htbp]
\centering
\subfigure[Schematic of target configuration.]{\includegraphics[width=0.45\textwidth]{Pore3schematic.pdf}}\qquad
\subfigure[Splitting probabilities and dynamic fluxes.]{\includegraphics[width=0.45\textwidth]{splittingfig.pdf}}
\caption{Left: A schematic of the three target configuration from Sec.~\ref{sec:resultsPlanar}. Right: The fractional fluxes $q_k(t) = \mathcal{J}_k(t)/\sum_{j=1}^N \mathcal{J}_j(t)$ into each target together with the steady state splitting probabilities $\phi_k$ obtained from \eqref{eqn:splitsol}. \label{fig:splits}}
\end{figure}
To demonstrate this theory, we compare these static splitting probabilities with the time-dependent fractional fluxes $q_k(t) = \mathcal{J}_k(t)/\sum_{j=1}^N \mathcal{J}_j(t)$ into each target (see Fig.~\ref{fig:splits}) obtained from the hybrid-asymptotic method (Sec.~\ref{sec:asy}). We draw attention to two important conclusions from this example. First, the dynamic fluxes converge to the static splitting probabilities on a very long timescale. For many physically or biologically relevant timescales, this questions the usefulness of using splitting probabilities for inference \cite{Miles2020} purposes. Second, the ordering of the relative fluxes into each target changes over the displayed time interval. Specifically, at short times ($t\lessapprox10^1$), target 1 captures almost all the flux and indeed is the most significant absorber over the entire timeline. Target 1 is the smallest target but closest to the initial location demonstrating that this distance is a significant indicator of capture ability. Later in the timeline, we see that targets 2 and 3 interchange their prominence in capture fluxes (around $t\approx 10^4$) implying that proximity to the initial location promotes faster capture at short times while at larger times, target size can play a more significant role. Importantly, none of these subtleties are apparent from the static splitting probabilities highlighting the necessity of obtaining full time dependent statistics.

\section{Discussion}
In this work we have demonstrated several methodologies for obtaining arrival time statistics of diffusing particles to complex sets of absorbing targets in planar regions. The Laplace transform approach seeks to solve an equation of modified Helmholtz type by either an asymptotic expansion for well separated targets or a boundary integral representation. In both cases the target geometries can be very general, but for the boundary integral method, it is presently limited to purely absorbing targets. The inverse transform is obtained by quadrature of the Bromwich integral. To complement these approaches, we developed a particle based kinetic Monte-Carlo (KMC) method that can resolve the arrival distribution for very general configurations of targets and boundary sets. These methods are rapid, accurate and easy to implement. The hybrid asymptotic method is particularly suited to the scenario of well-separated targets while the KMC method is applicable to general geometric scenarios and varied boundary conditions.

\begin{figure}
\centering
\includegraphics[width = 0.45\textwidth]{crazypath.pdf}
\caption{A tricky first passage time distribution to evaluate using either KMC or hybrid-asymptotic methods.\label{fig:crazypath}}
\end{figure}

There are a few conclusions that emerge from our study. Homogenization is a powerful technique that accurately reproduces the first passage time distributions of complex target sets by replacing them with a single circular target of appropriately chosen radius. However, homogenization brings limitations with it, particularly as it coarse grains the spatial distribution of arrivals across the targets. The relative fraction of particles that arrive across a distribution of targets has directional information that can be used to infer source location \cite{Miles2020}. Additionally, the dynamics available from the full distribution of arrival statistics reveals the limitations of using static information, (e.g.~splitting probabilities) which are only representative of very long time behavior. Our example of dynamic splitting probabilities (Sec.~\ref{sec:resultsSplitting}) suggests that there are several timescales over which arrival statistics can be important and that the relevant physical or biological timescale must be considered when inferring source location from arrival information.

While the methods developed here are quite general, there are certain scenarios where they fail and new approaches are needed. For example, in the scenario where there are numerous targets but only a few are reactive (Fig.~\ref{fig:crazypath}), particles must navigate a torturous route through inert targets to reach the destination. The combination of targets with either purely Neumann or Dirichlet boundary conditions hampers all the approaches developed here. This scenario is particularly challenging for the KMC method due to the fact that on a reflecting target surface, the particle will tend to perform a surface diffusion characterized by many small jumps. The very long wait time for a large jump necessary to leave the target vicinity means that the convergence rate is greatly reduced. To address these shortcoming, we plan extensions to the boundary integral formulation.

redpajama_set_name: RedPajamaArXiv


text:
You Belong: A Call for Connection (Paperback)
By Sebene Selassie
"A POWERFUL WORK OF SPIRITUALITY AND ANTI-RACISM"—Publishers Weekly
"IF YOU READ ONE BOOK IN 2020, MAKE IT THIS ONE."—Tricycle
From much-admired meditation expert Sebene Selassie, You Belong is a call to action, exploring our tangled relationship with belonging, connection, and each other
You are not separate. You never were. You never will be.
We are not separate from each other. But we don't always believe it, and we certainly don't always practice it. In fact, we often practice the opposite—disconnection and domination. From unconscious bias to "cancel culture," denial of our inherent interconnection limits our own freedom.
In You Belong, much-admired meditation expert Sebene Selassie reveals that accepting our belonging is the key to facing the many challenges currently impacting our world. Using ancient philosophy, multidisciplinary research, exquisite storytelling, and razor-sharp wit, Selassie leads us in an exploration of all the ways we separate (and thus suffer) and offers a map back to belonging.
To belong is to experience joy in any moment: to feel pleasure, dance in public, accept death, forgive what seems unforgivable, and extend kindness to yourself and others. To belong is also to acknowledge injustice, reckon with history, and face our own shadows. Full of practical advice and profound revelations, You Belong makes a winning case for resisting the forces that demand separation and reclaiming the connection—and belonging—that have been ours all along.
Sebene Selassie is a teacher & author who guides people to remember and trust their belonging. Born in Ethiopia and raised in Washington DC, she is trained as a meditation teacher, an integral coach, and as a practitioner of Indigenous Focusing Oriented Therapy (IFOT). She offers courses, workshops, and retreats online and in person and is one of the most popular teachers on the Ten Percent Happier meditation app. She lives with her partner Frederic on unceded Lenape territory in Brooklyn, NY.
"If you read one book in 2020, make it this one."
— Tricycle
"With grace, intelligence, and palpable compassion, Selassie weaves discussions of implicit bias, intergenerational trauma, and white supremacy throughout her coverage of more traditional spiritual topics, such as ego and conscious awareness. . . . This is a powerful work of spirituality and anti-racism."
— Publishers Weekly
"Sebene Selassie is one of the most highly respected and beloved newer voices in the mindfulness world. Her vision of belonging is one where we all find a home somewhere...in this body, this mind, with one another, on this planet. With depth and understanding (plus flashes of really funny humor) in this book Sebene skillfully traces a path we would all do well to follow."
— Sharon Salzberg, Bestselling author of Real Change
"This powerful and timely book shines a clear light on the root of our personal and societal suffering, the misperception that we don't belong. Sebene is a wise, brave and trustworthy spiritual guide—the teachings, stories and reflections in You Belong reveal the shared field of awareness and love that is our true home."
— Tara Brach, Author of Radical Compassion
"In You Belong, Sebene Selassie guides us in exploring all the ways we feel we don't belong. The good news is we do belong – to everyone and everything. That may be the bad news too. We belong to the feelings of loneliness, anxiety, and separation as much as to the love and connection. We belong to each other (even those we don't like or agree with), and our joy and freedom depend on recognizing this truth. With wisdom and humor, this book shows us that belonging has always been possible and it's the only way through our challenging times."
— Iyanla Vanzant, Bestselling author of Acts of Faith
Publisher: HarperOne
Publication Date: January 25th, 2022
Self-help / Motivational & Inspirational
Body, Mind & Spirit / Mindfulness & Meditation
Philosophy / Buddhist
Kobo eBook (August 25th, 2020): $12.99
Hardcover (August 25th, 2020): $27.99
Compact Disc (August 25th, 2020): $34.99
MP3 CD (August 25th, 2020): $39.99

redpajama_set_name: RedPajamaCommonCrawl


text:
More Light From Metamaterials
SemiconductorsNewsTypeTopic
Hyperbolic metasurface gives new control over lightwaves
Neil Savage
Image: Harvard University
nanotechnologydevicesHarvardquantum computinghyperbolic metasurfacemetamaterialsinvisibility
A material that manipulates light in unusual ways could lead to a whole variety of exotic devices, including microscopes capable of seeing inside cells, optical circuits for quantum computers, and invisibility cloaks.
The material in question is a hyperbolic metasurface, a two-dimensional type of metamaterial with a negative index of refraction that bends light in directions it would not normally travel, sending it along a hyperbola rather than an ellipse. One difficulty with metamaterials is that they often contain metals that absorb photons, limiting the distance light can travel to a few hundred nanometers. Diffraction, in which light bends around objects in its path, can also distort a lightwave in a metamaterial, limiting the reach of a light-based signal.
Hyperbolic metasurfaces could steer light according to its wavelength or by a magnetic field
This new metasurface overcomes the distance problem by sending light along the surface of a metal grating rather than through it, avoiding absorption. Harvard chemistry and physics professor Hongkun Park and his team describe the metasurface in the current issue of Nature.
They started by growing a single crystal of silver on top of a piece of silicon; then they used plasma to etch a grating into the silver. When they shone a laser onto the grating, surface plasmon polaritons—oscillations of electron density—formed at the interface of the silver and surrounding air, carrying the light along. The single crystal meant the surface was so smooth that there were no defects to absorb the light. The team's new technique increased the propagation distance of the light by two orders or magnitude, Park says. The light was also free of diffraction, which means it could be used to image objects much smaller than its wavelength.
Alexander High, a postdoctoral research in Park's group, says the metasurface is relatively easy to build. "There's no single aspect of it that is incredibly challenging or difficult to implement," he says.
The hyperbolic metasurface gives new levels of control over the propagation of light. The size of the grating determines which wavelengths are negatively refracted, so the metasurface can be used to route light depending on its wavelength. It also makes possible control of light through a property of quantum mechanics known as spin. Such spin control would allow the direction of light to be switched by applying a magnetic field. Finally, it provides a way of delivering light to other tiny optical elements, such as quantum dots.
It even gives fine-grained enough control that it would be possible to build single-photon transitors. Such devices could be elements of a quantum computer. What's more, the new metasurface is a step further on the road to invisibility cloaks; most of the demonstrations of optical cloaking have been at infrared or microwave wavelengths, whereas this would bring it into the visible spectrum. "There are lots of different possibilities," Park says.
The InstituteTopicEnergyArticleType
IEEE Discusses 6 Simple Solutions to Climate Change at COP27
EnergyTopicNewsType
Smaller, Cheaper Flow Batteries Throw Out Decades-Old Designs
Wearable Ultrasound Patch Images the Heart in Real Time
ComputingTopicMagazineTypeFeatureSpecial ReportsJanuary 2023
An IBM Quantum Computer Will Soon Pass the 1,000-Qubit Mark
ComputingTopicNewsType
Engineers Push Probabilistic Computers Closer to Reality
Cosmic Rays Are Quantum Computers' Kryptonite
SemiconductorsTopicMagazineTypeFeatureComputingSpecial ReportsJanuary 2023
Avicena's blue microLEDs are the dark horse in a race with Ayar Labs' laser-based system
Samuel K. Moore
Avicena's microLED chiplets could one day link all the CPUs in a computer cluster together.
Avicena
If a CPU in Seoul sends a byte of data to a processor in Prague, the information covers most of the distance as light, zipping along with no resistance. But put both those processors on the same motherboard, and they'll need to communicate over energy-sapping copper, which slow the communication speeds possible within computers. Two Silicon Valley startups, Avicena and Ayar Labs, are doing something about that longstanding limit. If they succeed in their attempts to finally bring optical fiber all the way to the processor, it might not just accelerate computing—it might also remake it.
Both companies are developing fiber-connected chiplets, small chips meant to share a high-bandwidth connection with CPUs and other data-hungry silicon in a shared package. They are each ramping up production in 2023, though it may be a couple of years before we see a computer on the market with either product.
Ayar Labs, has succeeded at drastically miniaturizing and reducing the power consumption of the kinds of silicon-photonics components used today to sling bits around data centers through optical-fiber cables. That equipment encodes data onto multiple wavelengths of light from an infrared laser and sends the light through a fiber.
Avicena's chiplet couldn't be more different: Instead of infrared laser light, it uses ordinary light from a tiny display made of blue microLEDs. And instead of multiplexing all the optical data so it can travel down a single fiber, Avicena's hardware sends data in parallel through the separate pathways in a specialized optical cable.
Ayar has the weight of history on its side, offering customers a technology similar to what they already use to send data over longer distances. But Avicena, the dark horse in this race, benefits from ongoing advances in the microdisplay industry, which is predicted to grow 80 percent per year and reach US $123 billion by 2030, fueled by a future full of virtual-reality gear and even augmented-reality contact lenses.
"Those companies are two ends of the spectrum in terms of the risk and innovation," says Vladimir Kozlov, founder and CEO of LightCounting, a telecommunications analysis firm.
MicroLEDs vs. Infrared Lasers
Avicena's silicon chiplet, LightBundle, consists of an array of gallium-nitride microLEDs, an equal-size array of photodetectors, and some I/O circuitry to support communication with the processor it feeds with data. Twin 0.5-millimeter-diameter optical cables link the microLED array on one chiplet to the photodetectors on another and vice versa. These cables—similar to the imaging cables in some endoscopes—contain a bundle of fiber cores that line up with the on-chip arrays, giving each microLED its own light path.
Besides the existence of this type of cable, Avicena needed two other things to come together, explains Bardia Pezeshki, the company's CEO. "The first one, which I think was the most surprising to anyone in the industry, is that LEDs could be run at 10 gigabits per second," he says. "That is stunning" considering that the state of the art for visible-light communication systems just five years ago was in the hundreds of megahertz. But in 2021, Avicena researchers revealed a version of the microLED they dubbed cavity-reinforced optical micro-emitters, or CROMEs. The devices are microLEDs that have been optimized for switching speed by minimizing capacitance and sacrificing some efficiency at converting electrons to light.
Gallium nitride isn't something that's typically integrated on silicon chips for computing, but thanks to advances in the microLED-display industry, doing so is essentially a solved problem. Seeking bright emissive displays for AR/VR and other things, tech giants such as Apple, Google, and Meta have spent years coming up with ways to transfer already-constructed micrometer-scale LEDs to precise spots on silicon and other surfaces. Now "it's done by the millions every day," says Pezeshki. Avicena itself recently purchased the fab where it developed the CROMEs from its Silicon Valley neighbor Nanosys.
Computer makers will want solutions that will not just help in the next two to three years but will give reliable improvements for decades.
The second component was the photodetector. Silicon isn't good at absorbing infrared light, so the designers of silicon-photonics systems typically compensate by making photodetectors and other components relatively large. But because silicon readily soaks up blue light, photodectors for Avicena's system need only be a few tenths of a micrometer deep, allowing them to be easily integrated in the chiplet under the imaging-fiber array. Pezeshki credits Stanford's David A.B. Miller with proving, more than a decade ago, that blue-light-detecting CMOS photodetectors were fast enough to do the job.
The combination of imaging fiber, blue microLEDs, and silicon photodetectors leads to a system that in prototypes transmits "many" terabits per second, says Pezeshki. Equally important as the data rate is the low energy needed to move a bit. "If you look at silicon-photonics target values, they are a few picojoules per bit, and these are from companies that are way ahead of us" in terms of commercialization, says Pezeshki. "We've already beaten those records." In a demo, the system moved data using about half a picojoule per bit. The startup's first product, expected in 2023, will not reach all the way to the processor but will aim to connect servers within a data-center rack. A chiplet for chip-to-chip optical links will follow "right on its heels," says Pezeshki.
But there are limits to the ability of microLEDs to move data. Because the LED light is incoherent, it suffers from dispersion effects that restrict it to about 10 meters. Lasers, in contrast, are naturally good at going the distance; Ayar's TeraPHY chiplets have a reach of up to 2 kilometers, potentially disrupting the architecture of supercomputers and data centers even more than Avicena's tech could. They could let computer makers completely rethink their architectures, allowing them to construct "essentially a single computer chip, but building it at rack scale," says Ayar CEO Charlie Wuischpard. The company is ramping up production with its partner GlobalFoundries and is building prototypes with partners in 2023, though these are not likely to be made public, he says.
Kozlov says to expect many more competitors to emerge. Computer makers will want solutions that will "not just help in the next two to three years but will give reliable improvements for decades." After all, the copper connections they are seeking to replace are still improving, too.
This article appears in the January 2023 print issue as "A Dark (Blue) Horse Emerges to Speed Up Computing."
HPE Invents First Memristor Laser ›
Plasmonics: A New Way to Link Processors With Light ›
Linking Chips With Light ›
Processor With Photonic Interconnects Built ›
Device Requirements for Optical Interconnects to Silicon Chips ... ›
CMOS-compatible multiple-wavelength oscillator for on-chip optical ... ›
Optical interconnect - Wikipedia ›