Optimization of the Mean First Passage Time in Near-Disk and Elliptical Domains in 2-D with Small Absorbing Traps

The concept of first passage time arises in various applications in biology, biochemistry, ecology, physics, and biophysics (see [6], [7], [20], [15] [23], [21], and the references therein). Narrow escape or capture problems are first passage time problems that characterize the expected time it takes for a Brownian “particle” to reach some absorbing set of small measure. These problems are of singular perturbation type as they involve two spatial scales: the ${\mathcal{O}}(1)$ spatial scale of the confining domain and the ${\mathcal{O}}(\varepsilon)$ asymptotically small scale of the absorbing set. Narrow escape and capture problems arise in various applications, including estimating the time it takes for a receptor to hit a certain target binding site, the time it takes for a diffusing surface-bound molecules to reach a localized signaling region on the cell membrane, or the time it takes for a predator to locate its prey, among others (cf. [1], [2], [4], [3], [9], [16], [24], [19], [15]). A comprehensive overview of the applications of narrow escape and capture problems in cellular biology is given in [8].

In this paper, we consider a narrow capture problem that involves determining the MFPT for a Brownian particle, confined in a bounded two-dimensional domain, to reach one of $m$ small stationary circular absorbing traps located inside the domain. The average MFPT for this diffusion process is the expected time for capture given a uniform distribution of starting points for the random walk. In the limit of small trap radius, this narrow capture problem can be analyzed by techniques in strong localized perturbation theory (cf. [26], [27]). For a disk-shaped domain spatial configurations of small absorbing traps that minimize the average MFPT domain were identified in [12]. However, the problem of identifying optimal trap configurations in other geometries is largely open. In this direction, the specific goal of this paper is to develop and implement a hybrid asymptotic-numerical theory to identify optimal trap configurations in near-disk domains and in the ellipse.

In § 2, we use a perturbation approach to derive a two-term approximation for the average MFPT in a class of near-disk domains in terms of a boundary deformation parameter $\sigma\ll 1$. In our analysis, we allow for a smooth, but otherwise arbitrary, star-shaped perturbation of the unit disk that preserves the domain area. At each order in $\sigma$, an approximate solution is derived for the MFPT that is accurate to all orders in $\nu\equiv{-1/\log\varepsilon}$, where $\varepsilon\ll 1$ is the common radius of the $m$ circular absorbing traps contained in the domain. To leading-order in $\sigma$, this small-trap singular perturbation analysis is formulated in the unit disk and leads to a linear algebraic system for the leading-order average MFPT involving the Neumann Green’s matrix. At order ${\mathcal{O}}(\sigma)$, a further linear algebraic system that sums all logarithmic terms in $\nu$ is derived that involves the Neumann Green’s matrix and certain weighted integrals of the boundary profile characterizing the domain perturbation. In § 3, we show how to numerically implement this asymptotic theory by using the analytical expression for the Neumann Green’s function for the unit disk together with the trapezoidal rule to compute certain weighted integrals of the boundary profile with high precision. From this numerical implementation of our asymptotic theory, and combined with either a simple gradient descent procedure or a particle swarming approach [11], we can numerically identify optimal trap configurations that minimize the average MFPT in near-disk domains. In § 3.1, we illustrate our hybrid asymptotic-numerical framework by determining some optimal trap configurations in various specific near-disk domains.

For a general 2-D domain containing small absorbing traps, a singular perturbation analysis in the limit of small trap radii, related to that in [15], [4], [12], and [26], shows that the average MFPT is closely approximated by the solution to a linear algebraic system involving the Neumann Green’s matrix. The challenge in implementing this analytical theory is that, for an arbitrary 2-D domain, a full PDE numerical solution of the Neumann Green’s function and its regular part is typically required to calculate this matrix. However, for an elliptical domain, in (4.36) and (4.37) below, we provide a new explicit representation of this Neumann Green’s function and its regular part. These explicit formulae allow for a rapid numerical evaluation of the Neumann Green’s interaction matrix for a given spatial distribution of the centers of the circular traps in the ellipse. The linear algebraic system determining the average MFPT is then coupled to a gradient descent numerical procedure in order to readily identify optimal trap configurations that minimize the average MFPT in an ellipse. Although, a similar formula for the Neumann Green’s function has been derived previously for a rectangular domain (cf. [17], [18], [14]), and an explicit and simple formula exists for the disk [12], to our knowledge there has been no prior derivation of a rapidly converging infinite series representation for the Neumann Green’s function in an ellipse. The derivation of this Neumann Green’s function using elliptic cylindrical coordinates is deferred until § 5.

With this explicit approach to determine the Neumann Green’s matrix, in § 4 we develop a hybrid asymptotic-numerical framework to approximate optimal trap configurations that minimize the average MFPT in an ellipse of a fixed area. In § 4.1 we implement our hybrid method to investigate how the optimal trap patterns change as the aspect ratio of the ellipse is varied. The results from the hybrid theory for the ellipse are favorably compared with full PDE numerical results computed from a computationally intensive numerical procedure of using the closest point method [10] to compute the average MFPT and a particle swarming approach [11] to numerically identify the optimum trap configuration. As the ellipse becomes thinner, our hybrid theory shows that the optimal trap pattern for $m=2,\ldots,5$ identical traps becomes collinear along the semi-major axis of the ellipse. In the limit of a long and thin ellipse, in § 4.2 a thin-domain asymptotic analysis is formulated and implemented to accurately predict the optimal locations of collinear trap configurations and the corresponding optimal average MFPT.

In § 6, we show that the optimal trap configurations that minimize the average MFPT also correspond to trap patterns that maximize the coefficient of order ${\mathcal{O}}(\nu^{2})$ in the asymptotic expansion of the fundamental Neumann eigenvalue of the Laplacian in the perforated domain. This fundamental eigenvalue characterizes the rate of capture of the Brownian particle by the traps. Eigenvalue optimization problems for the fundamental Neumann eigenvalue in a domain with small absorbing traps have been studied in [12] for the unit disk. The results herein extend this previous analysis to the ellipse and to near-disk domains.

We derive an asymptotic approximation for the MFPT for a class of near-disk 2-D domains that are defined in polar coordinates by

(2.1)

$$\Omega_{\sigma}=\Big{\{}(r,\theta)\,\Big{|}\,0<r\leq 1+\sigma h(\theta)\,,\,\,0\leq\theta\leq 2\pi\Big{\}}\,,$$

where the boundary profile, $h(\theta)$, is assumed to be an ${\mathcal{O}}(1)$, $C^{\infty}$ smooth $2\pi$ periodic function with $\int_{0}^{2\pi}h(\theta)\,d\theta=0$. Observe that $\Omega_{\sigma}\to\Omega$ as $\sigma\to 0$, where $\Omega$ is the unit disk. Since $\int_{0}^{2\pi}h(\theta)\,d\theta=0$, the domain area $|\Omega_{\sigma}|$ for $\sigma\ll 1$ is $|\Omega_{\sigma}|=\pi+{\mathcal{O}}(\sigma^{2})$.

Inside the perturbed disk $\Omega_{\sigma}$, we assume that there are $m$ circular traps of a common radius $\varepsilon\ll 1$ that are centered at arbitrary locations $\mathbf{x}_{1},\ldots,\mathbf{x}_{m}$ with $|\mathbf{x}_{i}-\mathbf{x}_{j}|={\mathcal{O}}(1)$ and $\mbox{dist}(\partial\Omega_{\sigma},\mathbf{x}_{j})={\mathcal{O}}(1)$ as $\varepsilon\to 0$. The $j$-th trap, centered at some $\mathbf{x}_{j}\in\Omega_{\sigma}$, is labelled by $\Omega_{\varepsilon j}=\{\mathbf{x}:|\mathbf{x}-\mathbf{x}_{j}|\leq\varepsilon\}$. The near-disk domain with the union of the trap regions deleted is denoted by $\bar{\Omega}_{\sigma}$. In $\bar{\Omega}_{\sigma}$, it is well-known that the mean first passage time (MFPT) for a Brownian particle starting at a point $\mathbf{x}\in\bar{\Omega}_{\sigma}$ to be absorbed by one of the traps satisfies (cf. [20])

(2.2)

$$\begin{split}D\,\Delta u=-1\,,&\quad\mathbf{x}\in\bar{\Omega}_{\sigma}\,;\qquad\bar{\Omega}_{\sigma}\equiv\Omega_{\sigma}\setminus\cup_{j=1}^{m}\Omega_{\varepsilon j}\,,\\ \partial_{n}u=0\,,\quad\mathbf{x}\in\partial\Omega_{\sigma}\,;&\qquad u=0\,,\quad\mathbf{x}\in\partial\Omega_{\varepsilon j}\,,\quad j=1,\ldots,m\,.\end{split}$$

In terms of polar coordinates, the Neumann boundary condition in (2.2) becomes

(2.3)

$$\begin{split}u_{r}-&\frac{\sigma h_{\theta}}{(1+\sigma h)^{2}}u_{\theta}=0\quad\text{on}\quad r=1+\sigma h(\theta)\,.\end{split}$$

For an arbitrary arrangement $\{{\mathbf{x}_{1},\ldots,\mathbf{x}_{m}\}}$ of the centers of the traps, and for $\sigma\to 0$ and $\varepsilon\to 0$, we will derive a reduced problem consisting of two linear algebraic systems that provide an asymptotic approximation to the MFPT that has an error ${\mathcal{O}}(\sigma^{2},\varepsilon^{2})$. These linear algebraic systems involve the Neumann Green’s matrix and certain weighted integrals of the boundary profile $h(\theta)$.

To analyze (2.2), we use a regular perturbation series to approximate (2.2) for the near-disk domain to problems involving a unit disk. We expand the MFPT $u$ as

(2.4)

$$\begin{split}u=u_{0}+\sigma u_{1}+\ldots\,,\end{split}$$

and substitute it into (2.2) and (2.3). This yields the leading-order problem

(2.5)

$$\begin{split}D\,\Delta u_{0}=-1\,,&\quad\mathbf{x}\in\bar{\Omega}\,;\qquad\bar{\Omega}\equiv\Omega\setminus\cup_{j=1}^{m}\Omega_{\varepsilon j}\,,\\ \partial_{n}u_{0}=0\,,\quad\text{on}\quad r=1\,;&\qquad u_{0}=0\,,\quad\mathbf{x}\in\partial\Omega_{\varepsilon j}\,,\quad j=1,\ldots,m\,,\end{split}$$

together with the following problem for the next order correction $u_{1}$:

(2.6)

$$\begin{split}\Delta u_{1}=0\,,&\quad\mathbf{x}\in\bar{\Omega}\,;\qquad\partial_{r}u_{1}=-hu_{0rr}+h_{\theta}u_{0\theta}\,,\quad\text{on}\quad r=1\,;\\ \quad u_{1}&=0\,,\quad\mathbf{x}\in\partial\Omega_{\varepsilon j}\,,\quad j=1,\ldots,m\,.\end{split}$$

Observe that (2.5) and (2.6) are formulated on the unit disk and not on the perturbed disk. Assuming $\varepsilon^{2}\ll\sigma$, we use (2.4) and $|\Omega_{\sigma}|=|\Omega|+{\mathcal{O}}(\sigma^{2})$ to derive an expansion for the average MFPT, defined by $\overline{u}\equiv\frac{1}{|\bar{\Omega}_{\sigma}|}\int_{\bar{\Omega}_{\sigma}}u\,\text{d}\mathbf{x}$, in the form

(2.7)

$$\begin{split}\overline{u}=\frac{1}{|\Omega|}\int_{\Omega}u_{0}\,\text{d}\mathbf{x}+\sigma\left[\frac{1}{|\Omega|}\int_{\Omega}u_{1}\,\text{d}\mathbf{x}+\frac{1}{|\Omega|}\int_{0}^{2\pi}h(\theta)\,u_{0}|_{r=1}\,\text{d}\theta\right]+\mathcal{O}(\sigma^{2},\varepsilon^{2})\,,\end{split}$$

where $|\Omega|=\pi$ and $u_{0}|_{r=1}$ is the leading-order solution $u_{0}$ evaluated on $r=1$.

Since the asymptotic calculation of the leading-order solution $u_{0}$ by the method of matched asymptotic expansions in the limit $\varepsilon\to 0$ of small trap radius was done previously in [4] (see also [15] and [26]), we only briefly summarize the analysis here. In the inner region near the $j$-th trap, we define the inner variables $\mathbf{y}=\varepsilon^{-1}(\mathbf{x}-\mathbf{x}_{j})$ and $u_{0}(\mathbf{x})=v_{j}(\varepsilon\mathbf{y}+\mathbf{x}_{j})$ with $\rho=|\mathbf{y}|$, for $j=1,\ldots,m$. Upon writing (2.5) in terms of these inner variables, we have for $\varepsilon\to 0$ and for each $j=1,\ldots,m$ that

(2.8)

$$\begin{split}\Delta_{\rho}\,v_{j}&=0\,,\quad\rho>1\,;\qquad v_{j}=0\,,\quad\mbox{on}\,\,\,\rho=1\,,\end{split}$$

where $\Delta_{\rho}\equiv\partial_{\rho\rho}+\rho^{-1}\partial_{\rho}$. This admits the radially symmetric solution $v_{j}=A_{j}\log\rho$, where $A_{j}$ is an unknown constant. From an asymptotic matching of the inner and outer solutions we obtain the required singularity condition for the outer solution $u_{0}$ as $\mathbf{x}\to\mathbf{x}_{j}$ for $j=1,\ldots,m$. In this way, we obtain that $u_{0}$ satisfies

(2.9a)		$\displaystyle\Delta u_{0}=-{1/D}\,,\quad\mathbf{x}$	$\displaystyle\in\Omega\setminus\{{\mathbf{x}_{1},\ldots,\mathbf{x}_{m}\}}\,;\quad\partial_{r}u_{0}=0\,,\,\,\,\mathbf{x}\in\partial\Omega\,;$
(2.9b)		$\displaystyle u_{0}\sim A_{j}\log|\mathbf{x}-\mathbf{x}_{j}|$	$\displaystyle+A_{j}/\nu\quad\text{as}\quad\mathbf{x}\to\mathbf{x}_{j}\,,\qquad j=1,\ldots,m\,,$

where $\nu\equiv-1/\log\varepsilon$. In terms of the Delta distribution, (2.9) implies that

(2.10)

$$\Delta u_{0}=-\frac{1}{D}+2\pi\sum_{j=1}^{m}A_{j}\delta(\mathbf{x}-\mathbf{x}_{j})\,,\quad\mathbf{x}\in\Omega\,;\qquad\partial_{r}u_{0}=0\,,\,\,\,\mathbf{x}\in\partial\Omega\,.$$

By applying the divergence theorem to (2.10) over the unit disk we obtain that $\sum_{j=1}^{m}A_{j}={|\Omega|/(2\pi D)}$. The solution to (2.10) is represented as

(2.11)

$\displaystyle u_{0}=-2\pi\sum_{k=1}^{m}A_{k}G(\mathbf{x};\mathbf{x}_{k})+\overline{u}_{0}\,;\qquad\overline{u}_{0}=\frac{1}{|\Omega|}\int_{\Omega}u_{0}\,\text{d}\mathbf{x}\,,$

where $G(\mathbf{x};\mathbf{x}_{j})$ is the Neumann Green’s function for the unit disk, which satisfies

(2.12a)		$$\displaystyle\Delta G=\frac{1}{|\Omega|}-\delta(\mathbf{x}-\mathbf{x}_{j})\,,\quad\mathbf{x}\in\Omega\,;\quad\partial_{n}G=0\,,\,\,\,\mathbf{x}\in\partial\Omega\,;\quad\int_{\Omega}G\,\text{d}\mathbf{x}=0\,,$$
(2.12b)		$$\displaystyle G\sim-\frac{1}{2\pi}\log{|\mathbf{x}-\mathbf{x}_{j}|}+R_{j}+\nabla_{\mathbf{x}}R_{j}\cdot(\mathbf{x}-\mathbf{x}_{j})\quad\text{as}\quad\mathbf{x}\to\mathbf{x}_{j}\,.$$

Here, $R_{j}\equiv R(\mathbf{x}_{j})$ is the regular part of the Green’s function at $\mathbf{x}=\mathbf{x}_{j}$. Expanding (2.11) as $\mathbf{x}\to\mathbf{x}_{j}$, and using the singularity behaviour of $G(\mathbf{x};\mathbf{x}_{j})$ given in (2.12b), together with the far-field behavior (2.9b) for $u_{0}$, we obtain the matching conditon:

(2.13)

$$-2\pi A_{j}\,R_{j}-2\pi\sum_{i\neq j}^{m}A_{i}\,G(\mathbf{x}_{j};\mathbf{x}_{i})+\overline{u}_{0}\sim{A_{j}/\nu}\,,\qquad\mbox{for}\quad j=1,\ldots,m\,.$$

This yields a linear algebraic system for $\overline{u}_{0}$ and $\mathcal{A}\equiv(A_{1},\ldots,A_{m})^{T}$, given by

(2.14)

$\displaystyle(I+2\pi\nu\,\mathcal{G})\mathcal{A}=\nu\,\overline{u}_{0}\,\mathbf{e}\,,\qquad\mathbf{e}^{T}\mathcal{A}=\frac{|\Omega|}{2\pi D}\,.$

Here, $\mathbf{e}\equiv(1,\ldots,1)^{T}$, $\nu=-1/\log\varepsilon$, $I$ is the $m\times m$ identity matrix, and $\mathcal{G}$ is the symmetric Green’s matrix with matrix entries given by

(2.15)

$\displaystyle(\mathcal{G})_{jj}=R_{j}\,\,\,\text{for}\,\,\,i=j\quad\text{and}\quad(\mathcal{G})_{ij}=(\mathcal{G})_{ji}=G(\mathbf{x}_{i};\mathbf{x}_{j})\,\,\,\text{for}\,\,\,i\neq j\,.$

We left-multiply the equation for $\mathcal{A}$ in (2.14) by $\mathbf{e}^{T}$, which isolates $\overline{u}_{0}$. By using this expression in (2.14), and defining the matrix $E$ by $E={\mathbf{e}\mathbf{e}^{T}/m}$, we get

(2.16)

$$\Big{[}I+2\pi\nu(I-E)\mathcal{G}\Big{]}\mathcal{A}=\frac{|\Omega|}{2\pi Dm}\mathbf{e}\,,\quad\text{and}\quad\overline{u}_{0}=\frac{|\Omega|}{2\pi D\nu m}+\frac{2\pi}{m}\mathbf{e}^{T}\mathcal{G}\mathcal{A}\,.$$

Next, we study the $\mathcal{O}(\sigma)$ problem for $u_{1}$ given in (2.6). We construct an inner region near each of the traps by introducing the inner variables $\mathbf{y}=\varepsilon^{-1}(\mathbf{x}-\mathbf{x}_{j})$ and $u_{1}(\mathbf{x})=V_{j}(\varepsilon\mathbf{y}+\mathbf{x}_{j})$ with $\rho=|\mathbf{y}|$. From (2.6), this yields the same leading-order inner problem (2.8) with $v_{j}$ replaced by $V_{j}$. The radially symmetric solution is $V_{j}=B_{j}\log\rho$, where $B_{j}$ is a constant to be found. By matching this far-field behavior of the inner solution to the outer solution we obtain the singularity behavior for $u_{1}$ as $\mathbf{x}\to\mathbf{x}_{j}$ for $j=1,\ldots,m$. In this way, we find from (2.6) that $u_{1}$ satisfies

(2.18a)		$\displaystyle\Delta u_{1}$	$\displaystyle=0\,,\quad\mathbf{x}\in\Omega\setminus\{\mathbf{x}_{1},\ldots,\mathbf{x}_{m}\}\,;\quad\partial_{r}u_{1}=F(\theta)\,,\quad\text{on}\quad r=1;$
(2.18b)		$\displaystyle u_{1}$	$\displaystyle\sim B_{j}\log{|\mathbf{x}-\mathbf{x}_{j}|}+B_{j}/\nu\quad\text{as}\quad\mathbf{x}\to\mathbf{x}_{j}\quad j=1,\ldots,m\,,$
where $\nu=-1/\log\varepsilon$ and $F(\theta)$ is defined by
(2.18c)		$$F(\theta)\equiv-hu_{0rr}|_{r=1}+h_{\theta}u_{0\theta}|_{r=1}=\left(hu_{0\theta}\right)_{\theta}+\frac{h}{D}\,.$$
In deriving (2.18c) we used $u_{0rr}=-u_{0\theta\theta}+{1/D}$ at $r=1$, as obtained from (2.5).

Next, we introduce the Dirac distribution and write the problem (2.18) for $u_{1}$ as

(2.19)

$\displaystyle\Delta u_{1}=2\pi\sum_{i=1}^{m}B_{i}\,\,\delta(\mathbf{x}-\mathbf{x}_{i})\,,\quad\mathbf{x}\in\Omega\,;\qquad u_{1r}=F(\theta)\,,\quad\text{on}\quad r=1\,.$

Since $\int_{0}^{2\pi}F(\theta)\,\text{d}\theta=0$, the divergence theorem yields $\sum_{j=1}^{m}B_{j}=0$. We decompose

(2.20)

$$u_{1}=-2\pi\sum_{i=1}^{m}B_{i}G(\mathbf{x};\mathbf{x}_{i})+u_{1p}+\overline{u}_{1}\,,$$

where $\overline{u}_{1}$ is the unknown average of $u_{1}$ over the unit disk, and $G(\mathbf{x};\mathbf{x}_{i})$ is the Neumann Green’s function satisfying (2.12). Here, $u_{1p}$ is taken to be the unique solution to

(2.21)

$\displaystyle\Delta u_{1p}$

$\displaystyle=0,\quad\mathbf{x}\in\Omega;\quad\partial_{r}u_{1p}=F(\theta)\,\quad\text{on}\quad r=1;\quad\int_{\Omega}u_{1p}\,\text{d}\mathbf{x}=0\,.$

Next, by expanding (2.20) as $\mathbf{x}\to\mathbf{x}_{j}$, we use the singularity behaviour of $G(\mathbf{x};\mathbf{x}_{j})$ as given in (2.12b) to obtain the local behavior of $u_{1}$ as $\mathbf{x}\to\mathbf{x}_{j}$, for each $j=1,\ldots,m$. The asymptotic matching condition is that this behavior must agree with that given in (2.18b). In this way, we obtain a linear algebraic system for the constant $\overline{u}_{1}$ and the vector $\mathbf{B}=(B_{1},\ldots,B_{m})^{T}$, which is given in matrix form by

(2.22)

$$(I+2\pi\nu\mathcal{G})\mathbf{B}=\nu\overline{u}_{1}\mathbf{e}+\nu\mathbf{u}_{1p}\,,\qquad\mathbf{e}^{T}\mathbf{B}=0\,.$$

Here, $I$ is the identity, $\mathbf{e}=(1,\ldots,1)^{T}$, and $\mathbf{u}_{1p}=(u_{1p}(\mathbf{x}_{1}),\ldots,u_{1p}(\mathbf{x}_{m}))^{T}$. Next, we left multiply the equation for $\mathbf{B}$ by $\mathbf{e}^{T}$. This determines $\overline{u}_{1}$, which is then re-substituted into (2.22) to obtain the uncoupled problem

(2.23)

$$\Big{[}I+2\pi\nu(I-E)\mathcal{G}\Big{]}\mathbf{B}=\nu(I-E)\mathbf{u}_{1p}\,,\quad\text{and}\quad\overline{u}_{1}=-\frac{1}{m}\mathbf{e}^{T}\mathbf{u}_{1p}+\frac{2\pi}{m}\mathbf{e}^{T}\mathcal{G}\mathbf{B}\,,$$

where $E\equiv{\mathbf{e}\mathbf{e}^{T}/m}$. Since $\mathbf{e}^{T}(I-E)=0$, we observe from (2.23) that $\mathbf{e}^{T}\mathbf{B}=0$, as required. Equation (2.23) gives a linear system for the $\mathcal{O}(\sigma)$ average MFPT $\overline{u}_{1}$ in terms of the Neumann Green’s matrix $\mathcal{G}$, and the vector $\mathbf{u}_{1p}$.

To determine $u_{1p}(\mathbf{x}_{j})$, we use Green’s second identity on (2.21) and (2.12) to obtain a line integral over the boundary $\mathbf{x}\in\partial\Omega$ of the unit disk. Then, by using (2.18c) for $F(\theta)$, integrating by parts and using $2\pi$ periodicity we get

(2.24)

$$u_{1p}(\mathbf{x}_{j})=\int_{0}^{2\pi}G(\mathbf{x};\mathbf{x}_{j})F(\theta)\,d\theta=\int_{0}^{2\pi}G(\mathbf{x};\mathbf{x}_{j})\frac{h(\theta)}{D}\,d\theta-\int_{0}^{2\pi}h(\theta)u_{0\theta}\partial_{\theta}G(\mathbf{x};\mathbf{x}_{j})\,d\theta\,.$$

Then, by setting (2.11) for $u_{0}$ into (2.24), we obtain in terms of the $A_{k}$ of (2.16) that

(2.25a)		$$u_{1p}(\mathbf{x}_{j})=\frac{1}{D}\int_{0}^{2\pi}G(\mathbf{x};\mathbf{x}_{j})h(\theta)\,d\theta+2\pi\sum_{k=1}^{m}A_{k}J_{jk}\,.$$
Here, $J_{jk}$ is defined by the following boundary integral with $\mathbf{x}=(\cos(\theta),\sin(\theta))^{T}$:
(2.25b)		$$J_{jk}\equiv\int_{0}^{2\pi}h(\theta)\left(\partial_{\theta}G(\mathbf{x};\mathbf{x}_{j})\right)\left(\partial_{\theta}G(\mathbf{x};\mathbf{x}_{k})\right)\,d\theta\,.$$

¿From a numerical evaluation of the boundary integrals in (2.25), we can calculate $\mathbf{u}_{1p}=(u_{1p}(\mathbf{x}_{1}),\ldots,u_{1p}(\mathbf{x}_{m}))^{T}$, which specifies the right-hand side of the linear system (2.23) for $\mathbf{B}$. After determining $\mathbf{B}$, we obtain $\overline{u}_{1}$ from the second relation in (2.23). Finally, by substituting (2.11) for $u_{0}$ into (2.7), and recalling that $\int_{0}^{2\pi}h(\theta)\,d\theta=0$, we obtain a two-term expansion for the average MFPT given by

(2.26)

$$\overline{u}\sim\overline{u}_{0}+\sigma\left(\overline{u}_{1}-2\sum_{k=1}^{m}A_{k}\int_{0}^{2\pi}G(\mathbf{x};\mathbf{x}_{k})h(\theta)\,d\theta\right)\,.$$

Here, $\mathbf{x}\in\partial\Omega$ and $\overline{u}_{0}$ is determined from (2.16).

To numerically evaluate the boundary integrals in (2.25) and (2.26), we need explicit formulae for $G(\mathbf{x};\mathbf{x}_{j})$ and $\partial_{\theta}G(\mathbf{x};\mathbf{x}_{j})$ on the boundary of the unit disk where $\mathbf{x}=(\cos\theta,\sin\theta)^{T}$. For the unit disk, we obtain from equation (4.3) of [12] that

(3.27a)		$$\displaystyle G(\mathbf{x};\mathbf{x}_{j})=-\frac{1}{2\pi}\log|\mathbf{x}-\mathbf{x}_{j}|-\frac{1}{4\pi}\log\left(|\mathbf{x}|^{2}|\mathbf{x}_{j}|^{2}+1-2\mathbf{x}\cdot\mathbf{x}_{j}\right)+\frac{(|\mathbf{x}|^{2}+|\mathbf{x}_{j}|^{2})}{4\pi}-\frac{3}{8\pi},$$
(3.27b)		$$\displaystyle R(\mathbf{x}_{j};\mathbf{x}_{j})=-\frac{1}{2\pi}\log\left(1-|\mathbf{x}_{j}|^{2}\right)+\frac{|\mathbf{x}_{j}|^{2}}{2\pi}-\frac{3}{8\pi}\,.$$

For an arbitrary configuration $\{{\mathbf{x}_{1},\ldots,\mathbf{x}_{m}\}}$ of traps, these expressions can be used to evaluate the Neumann Green’s matrix $\mathcal{G}$ of (2.15) as needed in (2.16) and (2.23).

Next, by setting $\mathbf{x}=(\cos\theta,\sin\theta)^{T}$ we can evaluate $G(\mathbf{x};\mathbf{x}_{j})$ on $\partial\Omega$, and then calculate its tangential boundary derivative $\partial_{\theta}G(\mathbf{x};\mathbf{x}_{j})$. By using (3.27a), we obtain

(3.28a)		$\displaystyle G(\mathbf{x};\mathbf{x}_{j})$	$\displaystyle=-\frac{1}{2\pi}\log\left(1+r_{j}^{2}-2r_{j}\cos(\theta-\theta_{j})\right)+\frac{1}{4\pi}(1+r_{j}^{2})-\frac{3}{8\pi}\,,$
(3.28b)		$\displaystyle\partial_{\theta}G(\mathbf{x};\mathbf{x}_{j})$	$\displaystyle=-\frac{r_{j}}{\pi}\frac{\sin(\theta-\theta_{j})}{\left[r_{j}^{2}+1-2r_{j}\cos(\theta-\theta_{j})\right]}\,,$

where $r_{j}\equiv|\mathbf{x}_{j}|$ and $\mathbf{x}_{j}=r_{j}(\cos\theta_{j},\sin\theta_{j})^{T}$. Then, since $\int_{0}^{2\pi}h(\theta)\,d\theta=0$, we can write the two boundary integrals appearing in (2.25) and (2.26) explicitly as

(3.29a)		$$\displaystyle\int_{0}^{2\pi}G(\mathbf{x};\mathbf{x}_{j})h(\theta)\,d\theta=-\frac{1}{2\pi}\int_{0}^{2\pi}h(\theta)\log\left(1+r_{j}^{2}-2r_{j}\cos(\theta-\theta_{j})\right)\,d\theta\,,$$
(3.29b)		$$\displaystyle J_{jk}=\frac{r_{j}r_{k}}{\pi^{2}}\int_{0}^{2\pi}\frac{h(\theta)\sin(\theta-\theta_{j})\sin(\theta-\theta_{k})}{\left[r_{j}^{2}+1-2r_{j}\cos(\theta-\theta_{j})\right]\left[r_{k}^{2}+1-2r_{k}\cos(\theta-\theta_{k})\right]}\,d\theta\,.$$

Although for an arbitrary $h(\theta)$ the integrals in (3.29) cannot be evaluated in closed form, they can be computed to a high degree of accuracy with relatively few grid points using the trapezoidal rule since this quadrature rule is exponentially convergent for $C^{\infty}$ smooth periodic functions [25]. When $|x_{j}|<1$, the logarithmic singularities off of the axis of integration for $J_{jk}$ in (3.29) are mild and pose no particular problem. In this way, we can numerically calculate the two-term expansion (2.26) for the average MFPT with high precision.

Then, to determine the optimal trap configuration we can either use the particle swarming approach [11], or the ODE relaxation dynamics scheme

(3.30)

$$\frac{d\mathbf{z}}{dt}=-\nabla_{\mathbf{z}}\overline{u}\,,\qquad\mbox{where}\quad\mathbf{z}\equiv(x_{1},y_{1},\ldots,x_{m},y_{m})^{T}\,,$$

and $\overline{u}$ is given in (2.26). Starting from an admissible initial state $\mathbf{z}|_{t=0}$, where $\mathbf{x}_{j}=(x_{j},y_{j})\in\Omega_{0}$ at $t=0$ for $j=1,\ldots,m$, the gradient flow dynamics (3.30) converges to a local minimum of $\overline{u}$. Because of our high precision in calculating $\overline{u}$, a centered difference scheme with mesh spacing $10^{-4}$ was used to estimate the gradient in (3.30).

Next, we consider the trap optimization problem in an ellipse of arbitrary aspect ratio, but with fixed area $\pi$. Our analysis uses a new explicit analytical formula, as derived in § 5, for the Neumann Green’s function $G(\mathbf{x};\mathbf{x}_{0})$ and its regular part $R_{e}$ of (5.54).

For $m$ circular traps each of radius $\varepsilon$, the average MFPT $\overline{u}_{0}$ satisfies (see (2.16))

(4.32)

$$\overline{u}_{0}=\frac{|\Omega|}{2\pi D\nu m}+\frac{2\pi}{m}\mathbf{e}^{T}\mathcal{G}\mathcal{A}\,,\quad\mbox{where}\quad\Big{[}I+2\pi\nu(I-E)\mathcal{G}\Big{]}\mathcal{A}=\frac{|\Omega|}{2\pi Dm}\mathbf{e}\,.$$

Here $E\equiv{\mathbf{e}\mathbf{e}^{T}/m}$, $\mathbf{e}=(1,\ldots,1)^{T}$, $\nu\equiv{-1/\log\varepsilon}$, and the Green’s matrix $\mathcal{G}$ depends on the trap locations $\{{\mathbf{x}_{1},\ldots,\mathbf{x}_{m}\}}$. To determine optimal trap configurations that are minimizers of the average MFPT, given in (4.32), we use the ODE relaxation scheme

(4.33)

$$\frac{d\mathbf{z}}{dt}=-\nabla_{\mathbf{z}}\overline{u}_{0}\,,\qquad\mbox{where}\quad\mathbf{z}\equiv(x_{1},y_{1},\ldots,x_{m},y_{m})\,.$$

In our implementation of (4.33), the gradient was approximated using a centered difference scheme with mesh spacing $10^{-4}$. The results shown below for the optimal trap patterns are confirmed from using a particle swarm approach [11].

The derivation of the Neumann Green’s function and its regular part in § 5 is based on mapping the elliptical domain to a rectangular domain using

(4.34a)		$$x=f\cosh\xi\cos\eta\,,\quad y=f\sinh\xi\sin\eta\,,\qquad f=\sqrt{a^{2}-b^{2}}\,.$$
With these elliptic cylindrical coordinates, the ellipse is mapped to the rectangle $0\leq\xi\leq\xi_{b}$ and $0\leq\eta\leq 2\pi$, where $a=f\cosh\xi_{b}$ and $b=f\sinh\xi_{b}$, so that
(4.34b)		$$f=\sqrt{a^{2}-b^{2}}\,,\qquad\xi_{b}=\tanh^{-1}\left(\frac{b}{a}\right)=-\frac{1}{2}\log\beta\,,\qquad\beta\equiv\left(\frac{a-b}{a+b}\right)\,.$$

To determine $(\xi,\eta)$, given a pair $(x,y)$, we invert the transformation (4.34a) using

(4.35a)		$$\xi=\frac{1}{2}\log\left(1-2s+2\sqrt{s^{2}-s}\right)\,,\quad s\equiv\frac{-\mu-\sqrt{\mu^{2}+4f^{2}y^{2}}}{2f^{2}}\,,\quad\mu\equiv x^{2}+y^{2}-f^{2}\,.$$
To recover $\eta$, we define $\eta_{\star}\equiv\sin^{-1}(\sqrt{p})$ and use
(4.35b)		$$\eta=\begin{cases}\eta_{\star},&\text{if }x\geq 0\,,\,\,y\geq 0\\ \pi-\eta_{\star},&\text{if }x<0\,,\,\,y\geq 0\\ \pi+\eta_{\star},&\text{if }x\leq 0\,,\,\,y<0\\ 2\pi-\eta_{\star},&\text{if }x>0\,,\,\,y<0\end{cases}\,,\quad\mbox{where}\quad p\equiv\frac{-\mu+\sqrt{\mu^{2}+4f^{2}y^{2}}}{2f^{2}}\,.$$

As derived in § 5, the matrix entries in $\mathcal{G}$ are obtained from the explicit result

(4.36a)		$$\begin{split}G(\mathbf{x};\mathbf{x}_{0})&=\frac{1}{4|\Omega|}\left(|\mathbf{x}|^{2}+|\mathbf{x}_{0}|^{2}\right)-\frac{3}{16|\Omega|}(a^{2}+b^{2})-\frac{1}{4\pi}\log\beta-\frac{1}{2\pi}\xi_{>}\\ &\qquad-\frac{1}{2\pi}\sum_{n=0}^{\infty}\log\left(\displaystyle\prod_{j=1}^{8}|1-\beta^{2n}z_{j}|\right)\,,\quad\mbox{for}\quad\mathbf{x}\neq\mathbf{x}_{0}\,,\end{split}$$
where $|\Omega|=\pi ab$, $\xi_{>}\equiv\max(\xi,\xi_{0})$, and the complex constants $z_{1},\ldots,z_{8}$ are defined in terms of $(\xi,\eta)$, $(\xi_{0},\eta_{0})$ and $\xi_{b}$ by
(4.36b)		$$\begin{split}z_{1}&\equiv e^{-|\xi-\xi_{0}|+i(\eta-\eta_{0})}\,,\quad z_{2}\equiv e^{|\xi-\xi_{0}|-4\xi_{b}+i(\eta-\eta_{0})}\,,\quad z_{3}\equiv e^{-(\xi+\xi_{0})-2\xi_{b}+i(\eta-\eta_{0})}\,,\\ z_{4}&\equiv e^{\xi+\xi_{0}-2\xi_{b}+i(\eta-\eta_{0})}\,,\quad z_{5}\equiv e^{\xi+\xi_{0}-4\xi_{b}+i(\eta+\eta_{0})}\,,\quad z_{6}\equiv e^{-(\xi+\xi_{0})+i(\eta+\eta_{0})}\,,\\ z_{7}&\equiv e^{|\xi-\xi_{0}|-2\xi_{b}+i(\eta+\eta_{0})}\,,\quad z_{8}\equiv e^{-|\xi-\xi_{0}|-2\xi_{b}+i(\eta+\eta_{0})}\,.\end{split}$$

Observe that the Dirac point at $\mathbf{x}_{0}=(x_{0},y_{0})$ is mapped to $(\xi_{0},\eta_{0})$. The transformation (4.34) and its inverse (4.35), determines $G(\mathbf{x};\mathbf{x}_{0})$ explicitly in terms of $\mathbf{x}\in\Omega$.