2.1. Inner solutions

The inner solution near each Dirichlet patch $\Gamma _{\varepsilon _j}$ can be found by introducing the local coordinates ${\boldsymbol{y}} = \varepsilon _j^{-1} \textbf{Q}_j ({{\boldsymbol{x}}}-{{\boldsymbol{x}}}_j)$ , where $\textbf{Q}_j$ is an appropriate rotation matrix to restrict ${\boldsymbol{y}} = (y_1,y_2)$ to the upper-half-plane ${\mathbb{H}}_2 = {\mathbb{R}} \times {\mathbb{R}}_+$ (the matrix $\textbf{Q}_j$ plays no role since $\textbf{Q}_j^\dagger \textbf{Q}_j = {\textbf{I}}$ and $|{\boldsymbol{y}}| = \varepsilon _j^{-1} |{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_j|$ ). We look for an inner solution near the patch $\Gamma _{\varepsilon _j}$ in the form

(2.2) \begin{equation} S_k({{\boldsymbol{x}}}_j + \varepsilon _j \textbf{Q}_j^\dagger {\boldsymbol{y}}) = \delta _{j,k} + A_j g_\infty (\,{\boldsymbol{y}}), \quad j\in \lbrace {1,\ldots ,N\rbrace }, \end{equation}

where $A_j$ is an unknown constant, and $g_\infty (\,{\boldsymbol{y}})$ is the Green’s function satisfying the canonical BVP given by

(2.3a) \begin{align} \Delta g_\infty &= 0 \quad \textrm {in}\, {\mathbb{H}}_2, \\[-8pt] \nonumber \end{align}

(2.3b) \begin{align} \partial _n g_\infty &= 0 \quad \textrm {on}\, y_2 = 0, \, |y_1| \geq 1; \qquad g_\infty = 0 \quad \textrm {on}\, y_2 =0, \, |y_1| \lt 1 , \\[-8pt] \nonumber \end{align}

(2.3c) \begin{align} g_\infty &\sim \ln |{\boldsymbol{y}}| + {\mathcal O}(1) \quad \textrm {as}\, |{\boldsymbol{y}}| \to \infty , \\[11pt] \nonumber \end{align}

(the subscript $\infty$ highlights infinite reactivity of the perfect patch, see below). The exact solution of this classical problem is given in Appendix A for completeness. The analysis below will require only the knowledge of the asymptotic behaviour of $g_\infty (\,{\boldsymbol{y}})$ at infinity. We recall that the constant term in this behaviour,

(2.4) \begin{equation} g_\infty (\,{\boldsymbol{y}}) \sim \ln |{\boldsymbol{y}}| - \ln (d) + o(1) \quad \textrm {as}\, |{\boldsymbol{y}}| \to \infty , \end{equation}

is determined by the logarithmic capacity $d$ of the interval $({-}1,1)$ , which is simply $d={1/2}$ . Now putting ${\boldsymbol{y}} = \varepsilon _j^{-1} \textbf{Q}_j ({{\boldsymbol{x}}}-{{\boldsymbol{x}}}_j)$ , we get that the far-field behaviour of the inner solution is

(2.5) \begin{equation} S_k \sim \delta _{j,k} + A_j \big [\ln |{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_j| - \ln (\varepsilon _j d) + o(1)\big] \quad \textrm {as}\ {{\boldsymbol{x}}}\to {{\boldsymbol{x}}}_j. \end{equation}

Setting

(2.6) \begin{equation} \nu _j = - \frac {1}{\ln (\varepsilon _j d)} , \quad \mbox{with} \quad d=\frac {1}{2} , \end{equation}

we rewrite this far-field behaviour, for each $j \in \lbrace {1,\ldots ,N\rbrace }$ , as

(2.7) \begin{equation} S_k \sim \delta _{j,k} + A_j \big[\ln |{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_j| + 1/\nu _j + o(1)\big] \quad \textrm {as}\,\,{{\boldsymbol{x}}}\to {{\boldsymbol{x}}}_j. \end{equation}

2.2. Outer solution and matching conditions

Now, we consider the outer problem for $S_k({{\boldsymbol{x}}})$ given by

(2.8a) \begin{align} \Delta S_k & = 0 \quad \textrm {in}\, \Omega , \\[-5pt] \nonumber \end{align}

(2.8b) \begin{align} S_k & \sim \delta _{j,k} + A_j\, \bigl [\ln |{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_j| + 1/\nu _j \bigr ] \quad \textrm {as}\ \ {{\boldsymbol{x}}} \to {{\boldsymbol{x}}}_j\in \partial \Omega , \quad j\in \lbrace {1,\ldots ,N\rbrace }, \\[-5pt] \nonumber \end{align}

(2.8c) \begin{align} \partial _n S_k & = 0 \quad \textrm {on}\, \partial \Omega \backslash \{{{\boldsymbol{x}}}_1,\ldots ,{{\boldsymbol{x}}}_N\}. \\[11pt] \nonumber \end{align}

To find this solution, we introduce the surface Neumann Green’s function $G({{\boldsymbol{x}}},{\boldsymbol \xi })$ , which satisfies

(2.9a) \begin{align} \Delta G & = \frac {1}{|\Omega |} \quad \textrm {in}\, \Omega ; \quad \partial _n G = 0 \quad \textrm {on}\, \partial \Omega \backslash \{{\boldsymbol \xi }\}; \quad \int \limits _{\Omega } G({{\boldsymbol{x}}},{\boldsymbol \xi })\,d{{\boldsymbol{x}}} = 0,\\[-10pt]\nonumber \end{align}

(2.9b) \begin{align} G({{\boldsymbol{x}}},{\boldsymbol \xi }) & \sim - \frac {1}{\pi } \ln |{{\boldsymbol{x}}}-{\boldsymbol \xi }| + R({\boldsymbol \xi }) + o(1) \quad \textrm {as}\, {{\boldsymbol{x}}} \to {\boldsymbol \xi }\in \partial \Omega , \end{align}

where $|\Omega |$ is the area of $\Omega$ , and $R({\boldsymbol \xi })$ is the regular part of $G({{\boldsymbol{x}}},{\boldsymbol \xi })$ , defined by

(2.10) \begin{equation} R({\boldsymbol \xi }) = \lim \limits _{{{\boldsymbol{x}}}\to {\boldsymbol \xi }} \biggl (G({{\boldsymbol{x}}},{\boldsymbol \xi }) + \frac {1}{\pi } \ln |{{\boldsymbol{x}}}-{\boldsymbol \xi }|\biggr ). \end{equation}

The formulation in (2.9) is equivalent to imposing that $\partial _n G=\delta ({{\boldsymbol{x}}}-{\boldsymbol \xi })$ on $\partial \Omega$ where ${\boldsymbol \xi }\in \partial \Omega$ .

The divergence theorem applied to Eq. (2.8) yields

(2.11) \begin{equation} \sum \limits _{j=1}^N A_j = 0. \end{equation}

Under this condition, the solution to Eq. (2.8) can be written as the linear combination

(2.12) \begin{equation} S_k({{\boldsymbol{x}}}) = \chi _k - \sum \limits _{i=1}^N \pi A_i G({{\boldsymbol{x}}},{{\boldsymbol{x}}}_i) , \end{equation}

where $\chi _k$ is a constant to be found. Matching the inner and outer asymptotic expansions as ${{\boldsymbol{x}}}\to {{\boldsymbol{x}}}_j$ gives

(2.13) \begin{equation} A_j + \nu _j\left[A_j R_j + \sum \limits _{i=1 \atop i\ne j}^N G_{j,i} A_i\right] = \chi _k \nu _j - \nu _k \delta _{j,k}, \end{equation}

for $j\in \lbrace {1,\ldots ,N\rbrace }$ , where we have defined

(2.14) \begin{equation} G_{j,i} = \pi G({{\boldsymbol{x}}}_j,{{\boldsymbol{x}}}_i) \quad (i\ne j), \qquad R_j = \pi R({{\boldsymbol{x}}}_j). \end{equation}

Together with the compatibility condition (2.11), we obtain a system of $N+1$ linear equations that determine the unknown coefficients $A_1, \ldots ,A_N$ and $\chi _k$ .

2.3. Matrix reformulation and general solution

To proceed, we introduce the following vectors and matrices of sizes $N\times 1$ and $N\times N$ :

(2.15) \begin{align} {\textbf{e}} = \left (\begin{array}{c} 1 \\ 1 \\ \cdots \\ 1 \\ \end{array}\right ), \quad {\boldsymbol \nu } = \left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \nu _1 & 0 & \cdots & 0 \\ 0 & \nu _2 & \cdots & 0 \\ \cdots &\cdots &\cdots &\cdots \\ 0 & 0 &\cdots & \nu _N \\ \end{array}\right ), \quad {\textbf{G}} = \left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} R_1 & G_{1,2} & \cdots & G_{1,N} \\ G_{2,1} & R_2 & \cdots & G_{2,N} \\ \cdots &\cdots &\cdots &\cdots \\ G_{N,1} & G_{N,2} &\cdots & R_N \\ \end{array}\right )\, \quad {\textbf{e}}_k = \left (\begin{array}{c} 0 \\ 1 \\ \cdots \\ 0 \\ \end{array}\right ),\nonumber\\[4pt] \end{align}

where $1$ stands on the $k$ -th row of the vector ${\textbf{e}}_k$ .

In terms of this notation, the system (2.13) is written in matrix form as

(2.16) \begin{equation} {\textbf{A}} + {\boldsymbol \nu } {\textbf{G}} {\textbf{A}} = \chi _k {\boldsymbol \nu } {\textbf{e}} - {\boldsymbol \nu } {\textbf{e}}_k . \end{equation}

Applying ${\textbf{e}}^\dagger$ on the left, we isolate $\chi _k$ as

(2.17) \begin{equation} \chi _k = \frac {{\textbf{e}}^\dagger {\boldsymbol \nu } {\textbf{G}} {\textbf{A}}}{\bar {\nu }} + \frac {\nu _k}{\bar {\nu }} , \end{equation}

where we used ${\textbf{e}}^\dagger {\textbf{A}} = 0$ due to Eq. (2.11), and defined

(2.18) \begin{equation} \bar {\nu } = {\textbf{e}}^\dagger {\boldsymbol \nu } {\textbf{e}} = \sum \limits _{j=1}^N \nu _j . \end{equation}

Eliminating $\chi _k$ from Eq. (2.16), we get a matrix equation for $\textbf{A}$ given by

(2.19) \begin{equation} ({\textbf{I}} + {\boldsymbol \nu } {\textbf{G}}) {\textbf{A}} - {\boldsymbol \nu } {\textbf{e}} \, \frac {{\textbf{e}}^\dagger {\boldsymbol \nu } {\textbf{G}} {\textbf{A}}}{\bar {\nu }} = \frac {\nu _k}{\bar {\nu }} {\boldsymbol \nu } {\textbf{e}} - \nu _k {\textbf{e}}_k . \end{equation}

By introducing the matrices ${\textbf{M}}_0$ and $\textbf{E}$ by

(2.20) \begin{equation} {\textbf{M}}_0 = {\textbf{I}} + \biggl ({\textbf{I}} - \frac {{\boldsymbol \nu } {\textbf{E}}}{\bar {\nu }}\biggr ) {\boldsymbol \nu } {\textbf{G}}, \qquad {\textbf{E}} = {\textbf{e}} {\textbf{e}}^\dagger , \end{equation}

the solution to Eq. (2.19) is

(2.21) \begin{equation} {\textbf{A}} = \frac {\nu _k}{\bar {\nu }} {\textbf{M}}_0^{-1} \bigl ({\boldsymbol \nu } {\textbf{e}} - \bar {\nu } {\textbf{e}}_k \bigr ), \end{equation}

where ${\textbf{M}}_0^{-1}$ is the inverse of ${\textbf{M}}_0$ . In fact, in the small-target limit, all $\nu _j \ll 1$ so that ${\textbf{M}}_0$ is a small perturbation of the identity matrix $\textbf{I}$ and is thus invertible. Once the coefficients $A_i$ are found, the constant $\chi _k$ follows from Eq. (2.17). As a consequence, the splitting probability $S_k({{\boldsymbol{x}}})$ is fully determined via the representation (2.12). This is the main result of this section. The shape of the confining domain $\Omega$ and the arrangement of Dirichlet patches are captured by the matrix $\textbf{G}$ , whereas the sizes of patches are accounted for via the matrix $\boldsymbol \nu$ . When the surface Neumann Green’s function is known analytically, a numerical computation of the coefficients $A_i$ and $\chi _k$ is fast, at least if the number of targets is not too large. We emphasise that the solution $\textbf{A}$ in (2.21) to the linear system (2.19) has accounted for all logarithmic correction terms in the asymptotic expansion of the splitting probability. This technique for effectively summing what otherwise would be an infinite logarithmic expansion in powers of $\nu _k$ was developed in [Reference Ward, Henshaw and Keller78], and has been used in other contexts (see [Reference Coombs, Straube and Ward20, Reference Kurella, Tzou, Coombs and Ward54, Reference Pillay, Ward, Peirce and Kolokolnikov65]).

We further emphasise that Eq. (2.12) is only applicable in the outer region, i.e., when $|{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_j| \gg {\mathcal O}(\varepsilon _j)$ for all $j \in \lbrace {1,\ldots ,N\rbrace }$ . In turn, if the starting point ${\boldsymbol{x}}$ is too close to ${{\boldsymbol{x}}}_j$ , this asymptotic formula may give wrong values (e.g., negative or exceeding $1$ ). In practice, the outer solution can be capped by $0$ and $1$ to avoid such invalid values, i.e., one can use $\max \{0, \min \{1, S_k({{\boldsymbol{x}}})\}\}$ instead of $S_k({{\boldsymbol{x}}})$ . We remark that if an accurate approximation of the splitting probability is needed near the patch, one has to use the corresponding inner solution.

We also note that the constant $\chi _k$ can be interpreted as the volume-averaged splitting probability. In fact, if the starting point ${\boldsymbol{x}}$ is not fixed but uniformly distributed in $\Omega$ , the average over the starting point yields

(2.22) \begin{equation} \overline {S}_k = \frac {1}{|\Omega |} \int \limits _{\Omega } S_k({{\boldsymbol{x}}}) \, d{{\boldsymbol{x}}} = \chi _k , \end{equation}

where we used $\int _{\Omega } G({{\boldsymbol{x}}},{{\boldsymbol{x}}}_i)\, d{{\boldsymbol{x}}}=0$ .

2.4. Example of two patches

In the case of two targets ( $N = 2$ ), the matrix ${\textbf{M}}_0$ from Eq. (2.20) reads

(2.23) \begin{equation} {\textbf{M}}_0 = {\textbf{I}} + \gamma \left (\begin{array}{c@{\quad}c} R_1 - G_{1,2} & G_{1,2} - R_2 \\[4pt] G_{1,2} - R_1 & R_2 - G_{1,2} \\ \end{array}\right ), \end{equation}

where $\gamma = \nu _1 \nu _2/(\nu _1 + \nu _2)$ . The inverse of this matrix is

(2.24) \begin{equation} {\textbf{M}}_0^{-1} = \frac {1}{1 + \gamma (R_1 + R_2 - 2G_{1,2})} \left ( {\textbf{I}} + \gamma \left (\begin{array}{c@{\quad}c} R_2 - G_{1,2} & R_2 -G_{1,2} \\[4pt] R_1 -G_{1,2} & R_1 - G_{1,2} \\ \end{array}\right ) \right ) . \end{equation}

Substituting this expression into Eq. (2.21), we find for $k = 1$ that

(2.25) \begin{equation} -A_1 = A_2 = \biggl (\frac {1}{\nu _1} + \frac {1}{\nu _2} + \bigl [R_1 + R_2 - 2G_{1,2}\bigr ]\biggr )^{-1} . \end{equation}

Then, by using Eq. (2.17) we determine $\chi _1$ as

(2.26) \begin{equation} \chi _1 = \frac {\nu _1 - A_2 \nu _1 (R_1 - G_{1,2}) + A_2\nu _2 (R_2 - G_{1,2})}{\nu _1+\nu _2} , \end{equation}

which can be further simplified as

(2.27) \begin{equation} \chi _1 = \frac {1/\nu _2 + (R_2 - G_{1,2})}{1/\nu _1 + 1/\nu _2 + (R_1 + R_2 - 2G_{1,2})} . \end{equation}

For instance, if $\Omega$ is the unit disk, the surface Neumann Green’s function is well known [Reference Pillay, Ward, Peirce and Kolokolnikov65]:

(2.28) \begin{equation} G({{\boldsymbol{x}}},{\boldsymbol \xi }) = - \frac {1}{\pi } \ln |{{\boldsymbol{x}}}-{\boldsymbol \xi }| + \frac {|{{\boldsymbol{x}}}|^2}{4\pi } - \frac {1}{8\pi } , \quad R({\boldsymbol \xi }) = \frac {1}{8\pi } . \end{equation}

Substitution of these expressions into Eqs. (2.25, 2.27) yields

(2.29a) \begin{align} -A_1 = A_2 & = \left({-}\ln (\varepsilon _1\varepsilon _2) + 2\ln (2) + 2\ln |{{\boldsymbol{x}}}_1-{{\boldsymbol{x}}}_2|\right)^{-1}, \\[-8pt] \nonumber \end{align}

(2.29b) \begin{align} \chi _1 & = \frac {-\ln (\varepsilon _2/2) + \ln |{{\boldsymbol{x}}}_1-{{\boldsymbol{x}}}_2|} {-\ln (\varepsilon _1\varepsilon _2/4) + 2\ln |{{\boldsymbol{x}}}_1-{{\boldsymbol{x}}}_2|} , \\[8pt] \nonumber \end{align}

and we conclude that

(2.30) \begin{equation} S_1({{\boldsymbol{x}}}) = \chi _1 + A_2 \ln \biggl (\frac {|{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_2|}{|{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_1|} \biggr ). \end{equation}

We can easily check that $S_1({{\boldsymbol{x}}})$ approaches $0$ [resp., $1$ ] as $\varepsilon _1\to 0$ [resp., $\varepsilon _2 \to 0$ ], as expected.

In the special case of two identical targets, $\nu _1 = \nu _2 = \nu$ , one has $\chi _1 = 1/2$ and $A_2 = 1/(2/\nu + 2\ln |{{\boldsymbol{x}}}_1-{{\boldsymbol{x}}}_2|)$ so that

(2.31) \begin{equation} S_1({{\boldsymbol{x}}}) = \frac 12 \biggl [1 + \frac {\nu }{1 + \nu \ln |{{\boldsymbol{x}}}_1-{{\boldsymbol{x}}}_2|} \ln \biggl (\frac {|{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_2|}{|{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_1|}\biggr )\biggr ]. \end{equation}

We remark that if we were to expand the denominator of the second term into a Taylor series in powers of $\nu \ll 1$ up to ${\mathcal O}(\nu ^2)$ , we would recover the truncated approximation given in Eq. (98) from [Reference Chevalier, Bénichou, Meyer and Voituriez14]. However, as $\nu = -1/\ln (\varepsilon /2)$ is not necessarily small enough, our new result Eq. (2.31) that incorporates all logarithmic terms is preferable to using the previous truncated approximation from [Reference Chevalier, Bénichou, Meyer and Voituriez14].

Figure 2 illustrates the splitting probability $S_1({{\boldsymbol{x}}})$ from Eq. (2.30). As explained earlier, we plot the capped version of this quantity, $\max \{0, \min \{1, S_1({{\boldsymbol{x}}})\}\}$ , to avoid invalid values near two patches. Expectedly, $S_1({{\boldsymbol{x}}})$ increases as ${\boldsymbol{x}}$ gets closer to the first patch (red arc) and decreases as ${\boldsymbol{x}}$ gets closer to the second patch (blue arc). However, as the outer solution (2.30) is not applicable in the vicinity of these patches, we can observe some discrepancy, e.g., $S_1({{\boldsymbol{x}}})$ does not vanish on the second patch, as it should. To amend this discrepancy, we can use the inner solution when $|{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_j| \lesssim \varepsilon _j$ .

Figure 2. Splitting probability $S_1({{\boldsymbol{x}}})$ , given by Eq. (2.30), for the unit disk with two Dirichlet patches of length $2\varepsilon _1 = 0.2$ (red) and $2\varepsilon _2 = 0.4$ (blue). Note that $S_1({{\boldsymbol{x}}})$ was capped by $0$ and $1$ , i.e., we plotted $\max \{0, \min \{1, S_1({{\boldsymbol{x}}})\}\}$ .

2.5. Example of equally spaced identical patches on the boundary of the unit disk

If there are $N$ identical targets, one has $\nu _j = \nu$ so that ${\boldsymbol \nu } = \nu {\textbf{I}}$ , $\bar {\nu } = N \nu$ , and thus ${\textbf{M}}_0 = {\textbf{I}} + \nu ({\textbf{I}} - {\textbf{E}}/N) {\textbf{G}}$ from Eq. (2.20). Moreover, if the patches are equally spaced on the boundary of the unit circle, the matrix $\textbf{G}$ defined in Eq. (2.15) is circulant and symmetric. As a consequence, its eigenvectors can be written as

(2.32) \begin{equation} \textbf{q}_j = \frac {1}{\sqrt {N}} \bigl (\omega ^{-j}, \omega ^{-2j}, \ldots , \omega ^{-Nj}\bigr )^\dagger , \qquad j\in \lbrace {1,\ldots ,N\rbrace }, \end{equation}

where $\omega = e^{2\pi i/N}$ , and the transposition $\dagger$ now denotes the Hermitian conjugate. Upon taking the real and imaginary parts of $\textbf{q}_j$ , the resulting real-valued eigenvectors form an orthonormal basis in ${\mathbb{R}}^{N}$ since $\textbf{G}$ is symmetric. Let us denote by $\kappa _j$ the associated eigenvalues of $\textbf{G}$ :

(2.33) \begin{equation} \kappa _j = R_1 + \sum \limits _{m=1}^{N-1} \omega ^{mj} \, G_{1,1+m}, \qquad j\in \lbrace {1,\ldots ,N\rbrace }. \end{equation}

For any $j \in \lbrace {1,\ldots ,N-1\rbrace }$ , we get

(2.34) \begin{equation} {\textbf{M}}_0 \textbf{q}_j = \bigl ({\textbf{I}} + \nu [{\textbf{I}} - {\textbf{E}}/N]\bigr ) \kappa _j \textbf{q}_j = (1 + \nu \kappa _j) \textbf{q}_j , \end{equation}

because ${\textbf{E}} \textbf{q}_j = {\textbf{e}} {\textbf{e}}^\dagger \textbf{q}_j = 0$ for any $1 \leq j \leq N-1$ due to orthogonality of $\textbf{q}_j$ to $\textbf{q}_N = {\textbf{e}}/\sqrt {N}$ . As a consequence, each $\textbf{q}_j$ with $j\in \lbrace {1,\ldots ,N-1\rbrace }$ is also the eigenvector of ${\textbf{M}}_0$ , associated to the eigenvalue $1 + \nu \kappa _j$ . In addition, we have

(2.35) \begin{equation} {\textbf{M}}_0 \textbf{q}_N = \bigl ({\textbf{I}} + \nu [{\textbf{I}} - {\textbf{E}}/N]\bigr ) \kappa _N \textbf{q}_N = \textbf{q}_N , \end{equation}

since ${\textbf{E}} \textbf{q}_N /N = \textbf{q}_N$ . Therefore, $\textbf{q}_N$ is the eigenvector of ${\textbf{M}}_0$ associated with the eigenvalue $1$ . We use this spectral information to invert the matrix ${\textbf{M}}_0$ as

(2.36) \begin{equation} {\textbf{M}}_0^{-1} = \textbf{q}_N \textbf{q}_N^\dagger + \sum \limits _{j=1}^{N-1} \textbf{q}_j (1 + \nu \kappa _j)^{-1} \textbf{q}_j^\dagger . \end{equation}

Substituting this spectral representation into Eq. (2.21), we get

(2.37) \begin{equation} {\textbf{A}} = - \nu \sum \limits _{j=1}^{N-1} \textbf{q}_j (1 + \nu \kappa _j)^{-1} \textbf{q}_j^\dagger {\textbf{e}}_k, \end{equation}

where we used the orthogonality of the eigenvectors $\textbf{q}_j$ . Substitution of this expression into Eq. (2.17) yields $\chi _k = 1/N$ . This is consistent with the interpretation of $\chi _k$ as the volume-averaged splitting probability: when all patches are identical and equally spaced on the boundary of the unit disk, they are equivalent from the uniformly distributed starting point, so that $\overline {S}_k = 1/N$ from Eq. (2.22).

To complete this example, we will simplify Eq. (2.33) by using the explicit form (2.28) of the surface Neumann Green’s function for the unit disk. Since the centres of the patches ${{\boldsymbol{x}}}_j$ are equally spaced on the domain boundary, we have ${{\boldsymbol{x}}}_j = e^{2\pi i(\,j-1)/N} = \omega ^{\,j-1}$ for $j\in \lbrace {1,\ldots ,N\rbrace }$ . In this way, substituting Eq. (2.28) into Eq. (2.33), we get

(2.38) \begin{equation} \kappa _j = \frac {1}{8}\sum _{m=0}^{N-1}\omega ^{\,jm} - \sum \limits _{m=1}^{N-1} \omega ^{mj} \ln \vert 1 - \omega ^m\vert , \qquad j\in \lbrace {1,\ldots ,N\rbrace } , \end{equation}

where we interpret points as complex numbers and $|z|$ as the modulus of $z$ . For $j=N$ , for which $\omega ^{Nm}=1$ , we get

(2.39) \begin{equation} \kappa _N = \frac {N}{8} - \ln \left \vert \prod _{m=1}^{N-1}(1-\omega ^m)\right \vert = \frac {N}{8}-\ln {N} , \end{equation}

where the product in Eq. (2.39) was evaluated by using the roots of unity together with L’Hopital’s rule to get $\lim _{z\to 1} {(z^N-1)/(z-1)}=N=\prod _{m=1}^{N-1}(1-\omega ^m)$ . For $j\lt N$ , the first term in Eq. (2.38) vanishes and we obtain

(2.40) \begin{equation} \kappa _j = - \sum \limits _{m=1}^{N-1} \omega ^{mj} \ln |1 - \omega ^m|, \qquad j\in \lbrace {1,\ldots ,N-1\rbrace } , \end{equation}

which reduces after some simplifications to

(2.41) \begin{equation} \kappa _j = \ln {2} - \sum _{m=1}^{N-1}\cos \left (\frac {2\pi j m}{N}\right ) \ln \left [ \sin \left (\frac {\pi m}{N}\right ) \right ], \end{equation}

for $j\in \lbrace {1,\ldots ,N-1\rbrace }$ . The asymptotic behaviour of $\kappa _j$ for large $N$ is derived in Appendix B.

3.1. The mean first-reaction time

Although our asymptotic analysis has focused on calculating splitting probabilities it can be easily modified to calculate the mean first-reaction time (MFRT).

The dimensionless MFRT $u({{\boldsymbol{x}}})$ satisfies a Poisson equation with mixed Neumann-Robin boundary conditions:

(3.18a) \begin{align} \Delta u & = - 1 , \quad {{\boldsymbol{x}}} \in \Omega , \\[-12pt] \nonumber \end{align}

(3.18b) \begin{align} \partial _{n} u + q_i u & = 0, \quad {{\boldsymbol{x}}} \in \Gamma _{\varepsilon _i} , \quad i\in \lbrace {1,\ldots ,N\rbrace } , \\[-12pt] \nonumber \end{align}

(3.18c) \begin{align} \partial _{n} u & = 0 , \quad {{\boldsymbol{x}}} \in \partial \Omega _0 =\partial \Omega \backslash (\Gamma _{\varepsilon _1} \cup \cdots \cup \Gamma _{\varepsilon _N}). \\[8pt] \nonumber \end{align}

As previously, each reactive boundary patch $\Gamma _{\varepsilon _i}$ has length $2\varepsilon _i$ , reactivity parameter $q_i$ , and is centred at ${{\boldsymbol{x}}}_i\in \partial \Omega$ .

The matched asymptotic analysis of Eq. (3.18) in the small-target limit $\varepsilon _i\ll 1$ is very similar to that for analysing the splitting probability. The inner solution near the $j$ -th patch in terms of an unknown coefficient $A_j$ is

(3.19) \begin{equation} V_{j}(\,{\boldsymbol{y}})= u({{\boldsymbol{x}}}_j+\varepsilon _j \textbf{Q}_j^\dagger {\boldsymbol{y}}) = A_j g_{\mu _j}(\,{\boldsymbol{y}}), \qquad j\in \lbrace {1,\ldots ,N\rbrace }, \end{equation}

where $\mu _j\equiv \varepsilon _j q_j$ and $g_{\mu }(\,{\boldsymbol{y}})$ is the Robin Green’s satisfying Eq. (3.3) with far-field behaviour Eq. (3.6). Upon matching the far-field behaviour of $V_{j}(\,{\boldsymbol{y}})$ to the outer solution, we find that to within all logarithmic terms the outer solution satisfies

(3.20a) \begin{align} \Delta u & = -1 \quad \textrm {in}\, \Omega , \\[-10pt] \nonumber \end{align}

(3.20b) \begin{align} u & \sim A_j \bigl [\ln |{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_j| + 1/\nu _j \bigr ] \quad \textrm {as}\, {{\boldsymbol{x}}} \to {{\boldsymbol{x}}}_j\in \partial \Omega , \quad j\in \lbrace {1,\ldots ,N\rbrace }, \\[-10pt] \nonumber\end{align}

(3.20c) \begin{align} \partial _n u & = 0 \quad \textrm {on}\, \partial \Omega \backslash \{{{\boldsymbol{x}}}_1,\ldots ,{{\boldsymbol{x}}}_N\}, \\[8pt] \nonumber \end{align}

where $\nu _j$ is defined by (3.12). The solvability condition for (3.20) is that $\sum _{j=1}^{N} A_j={|\Omega |/\pi }$ . We then represent $u$ in terms of the surface Neumann Green’s function and the volume average $\overline {u}_0=|\Omega |^{-1}\int _{\Omega } u({{\boldsymbol{x}}})\, d{{\boldsymbol{x}}}$ as

(3.21) \begin{equation} u({{\boldsymbol{x}}}) = \overline {u}_0 - \pi \sum \limits _{i=1}^N A_i G({{\boldsymbol{x}}},{{\boldsymbol{x}}}_i) . \end{equation}

Imposing the singularity behaviour in Eq. (3.20b), we obtain an $(N+1)$ -dimensional linear algebraic system for $\overline {u}_0$ and $A_1,\ldots ,A_N$ given by

(3.22) \begin{equation} A_j + \nu _j\left[A_j R_j + \sum \limits _{i=1 \atop i\ne j}^N G_{j,i} A_i\right] = \overline {u}_0 \nu _j, \quad j\in \lbrace {1,\ldots ,N\rbrace }; \qquad \sum _{j=1}^{N} A_j = \frac {|\Omega |}{\pi } , \end{equation}

where $G_{j,i}$ and $R_i$ were defined in (2.14). In matrix form Eq. (3.22) is written for ${\textbf{A}}\equiv (A_1,\ldots ,A_N)^\dagger$ as

(3.23) \begin{equation} {\textbf{A}} + {\boldsymbol \nu } {\textbf{G}} {\textbf{A}} = \overline {u}_0 {\boldsymbol \nu } {\textbf{e}} , \quad {\textbf{e}}^{T} {\textbf{A}} = \frac {|\Omega |}{\pi } , \end{equation}

where $\textbf{e}$ , $\boldsymbol \nu$ and the Green’s matrix $\textbf{G}$ were defined in Eq. (2.15). By eliminating $\overline {u}_0$ in Eq. (3.23), we conclude that

(3.24) \begin{equation} \overline {u}_0 = \frac {|\Omega |}{\pi \bar {\nu }} + \frac {{\textbf{e}}^\dagger {\boldsymbol \nu } {\textbf{G}} {\textbf{A}}}{\bar {\nu }} , \end{equation}

where $\bar {\nu }\equiv \sum \limits _{j=1}^N \nu _j$ , while $\textbf{A}$ is the solution to

(3.25) \begin{equation} {\textbf{M}}_0 {\textbf{A}} = \frac {|\Omega |}{\pi \bar {\nu }} {\boldsymbol \nu } , \end{equation}

with ${\textbf{M}}_0$ and $\textbf{E}$ being defined in Eq. (2.20).

By using Eq. (3.12) for $\nu _j$ , which involves the local reactivity parameter $q_j$ on the patch, one can invert ${\textbf{M}}_0$ in Eq. (3.25) to get the coefficients $A_j$ . As a consequence, Eq. (3.24) gives access to the volume-averaged MFRT $\overline {u}_0$ , whereas Eq. (3.21) determines the MFRT $u({{\boldsymbol{x}}})$ for any well-separated spatial configuration of partially reactive patches. This result generalises that in [Reference Pillay, Ward, Peirce and Kolokolnikov65], where perfect reactivities ( $q_j=\infty$ ) were assumed.

4.1. Matched asymptotic analysis

As before, we look at the inner solution near each Steklov patch $\Gamma _{\varepsilon _j}$ . Upon comparing the Robin and Steklov conditions (3.1, 4.1b), we notice two differences: (i) $q_j$ is replaced by $-\sigma$ , and (ii) there is no inhomogeneous term $q_j \delta _{j,k}$ in the right-hand side. Apart from these two points, the BVPs for $S_k({{\boldsymbol{x}}})$ and $V({{\boldsymbol{x}}})$ are identical. As a consequence, we can immediately rewrite the asymptotic behaviour (3.11) for each $j\in \lbrace {1,\ldots ,N\rbrace }$ as

(4.2) \begin{equation} V \sim A_j \left[\ln |{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_j| + 1/\nu _j + {\mathcal{C}}({-}\sigma \varepsilon _j) + o(1) \right] \quad \textrm {as}\, {{\boldsymbol{x}}}\to {{\boldsymbol{x}}}_j, \end{equation}

where we now redefine $\nu _j$ as

(4.3) \begin{equation} \nu _j = - {1/\ln (\varepsilon _j)} , \end{equation}

and where $A_j$ is an unknown coefficient. Note that the constant term ${\mathcal{C}}({-}\sigma \varepsilon _j)$ is not incorporated into the new definition of $\nu _j$ , as we did earlier in the Robin case.

As before, the outer solution is represented as

(4.4) \begin{equation} V = \chi - \pi \sum \limits _{i=1}^N A_i G({{\boldsymbol{x}}},{{\boldsymbol{x}}}_i), \end{equation}

with an unknown constant $\chi$ . The divergence theorem still ensures the compatibility condition (2.11). Upon enforcing the singularity behaviour (4.2) for the solution in Eq. (4.4), we get

(4.5) \begin{equation} \chi - ({\textbf{G}} {\textbf{A}})_j = A_j (1/\nu _j + {\mathcal{C}}({-}\sigma \varepsilon _j)), \quad j\in \lbrace {1,\ldots ,N\rbrace } , \end{equation}

where we used the matrix notations introduced in Eq. (2.15) of Section 2.3. Multiplying this equation by $\nu _j$ and introducing the diagonal matrix $\textbf{C}$ formed by $\{{\mathcal{C}}({-}\sigma \varepsilon _1), \ldots , {\mathcal{C}}({-}\sigma \varepsilon _N)\}$ , we rewrite Eq. (4.5) as

(4.6) \begin{equation} \bigl ({\textbf{I}} + {\boldsymbol \nu } {\textbf{G}} + {\boldsymbol \nu } {\textbf{C}}\bigr ) {\textbf{A}} = \chi {\boldsymbol \nu } {\textbf{e}} . \end{equation}

Left-multiplying this equation by ${\textbf{e}}^\dagger$ , and using ${\textbf{e}}^\dagger {\textbf{A}} = 0$ , we isolate $\chi$ as

(4.7) \begin{equation} \chi = \frac {1}{\bar {\nu }} \bigl ({\textbf{e}}^\dagger {\boldsymbol \nu }{\textbf{G}} + {\textbf{e}}^\dagger {\boldsymbol \nu }{\textbf{C}}\bigr ){\textbf{A}} , \end{equation}

where $\bar {\nu }$ was defined by Eq. (2.18). Substituting this expression back into Eq. (4.6), we obtain that

(4.8) \begin{equation} \biggl ({\textbf{I}} + \biggl ({\textbf{I}} - \frac {{\boldsymbol \nu } {\textbf{E}}}{\bar {\nu }}\biggr ) {\boldsymbol \nu } \left ( {\textbf{G}} + {\textbf{C}}\right ) \biggr ) {\textbf{A}} = \boldsymbol{0} . \end{equation}

The necessary and sufficient condition for the existence of a nontrivial solution to this matrix equation is

(4.9) \begin{equation} {\mathrm{det}}\biggl ({\textbf{I}} + \biggl ({\textbf{I}} - \frac {{\boldsymbol \nu } {\textbf{E}}}{\bar {\nu }}\biggr ) {\boldsymbol \nu } ({\textbf{G}} + {\textbf{C}})\biggr ) = \boldsymbol{0} . \end{equation}

The matrices $\boldsymbol \nu$ and $\textbf{G}$ are determined by the sizes and arrangement of the Steklov patches, while the matrix $\textbf{C}$ is formed by $\{{\mathcal{C}}({-}\sigma \varepsilon _1),\ldots , {\mathcal{C}}({-}\sigma \varepsilon _N)\}$ , with the function ${\mathcal{C}}(\mu )$ given by Eq. (3.7). As a consequence, Eq. (4.9) determines the unknown parameter $\sigma$ . Moreover, the functional form (3.7) implies that there are infinitely many negative solutions, denoted as $-\sigma _j$ , which are actually small-target approximations of the Steklov eigenvalues.

Let us first provide qualitative insights on these solutions. If the matrix $\textbf{C}$ was fixed in the small-target limit, the condition $\nu _j\ll 1$ would imply the smallness of the second matrix term in Eq. (4.9) as compared to the identity matrix $\textbf{I}$ , thus ensuring the positivity of the determinant. To compensate the smallness of the matrix $\boldsymbol \nu$ , the matrix $\textbf{C}$ must therefore be large in the small-target limit. This is possible when at least one $\sigma \varepsilon _i$ is close to $- \mu _{2k}$ for some $k$ . In other words, one can expect that a solution $-\sigma _j$ of Eq. (4.9) is close to $\mu _{2k}/\varepsilon _i$ for some $i$ and $k$ . This intuitive picture suggests that an eigenvalue $\sigma _j$ of the mixed Steklov-Neumann problem with $N$ patches can be approximated by that on a single patch (say, $\Gamma _{\varepsilon _i}$ ), as if there were no other patches and the associated eigenfunction was localised on $\Gamma _{\varepsilon _i}$ . The situation is, however, more subtle in the vicinity of $\mu _0 = 0$ . In the analysis below, we focus on the asymptotic behaviour of the first $N$ eigenvalues $\sigma _j$ .

4.2. First $N$ eigenvalues

Let us assume that $- \mu = \sigma \varepsilon _j \ll 1$ for all $j = 1,\ldots ,N$ , so that we can apply the approximate relation (3.10) for the function ${\mathcal{C}}(\mu )$ . Introducing the diagonal matrix $\boldsymbol \eta$ formed by $\{ {\pi /(2\varepsilon _1)}, \ldots , {\pi /(2\varepsilon _N)}\}$ , we get ${\textbf{C}} \approx - {\boldsymbol \eta }/\sigma + C_1 {\textbf{I}}$ . Its substitution into Eq. (4.9) yields

(4.10) \begin{equation} {\mathrm{det}}\biggl ({\textbf{I}} + \biggl ({\textbf{I}} - \frac {{\boldsymbol \nu } {\textbf{E}}}{\bar {\nu }}\biggr ) {\boldsymbol \nu } \bigl [{\textbf{G}} - {\boldsymbol \eta }/\sigma + C_1 {\textbf{I}}\bigr ]\biggr ) = \boldsymbol{0} . \end{equation}

Upon defining ${\textbf{M}}_1$ and ${\textbf{B}}_1$ by

(4.11) \begin{equation} {\textbf{M}}_1 = {\textbf{I}} + \biggl ({\textbf{I}} - \frac {{\boldsymbol \nu } {\textbf{E}}}{\bar {\nu }}\biggr ) {\boldsymbol \nu } \bigl ({\textbf{G}} + C_1 {\textbf{I}}\bigr ) , \qquad {\textbf{B}}_1 = \biggl ({\textbf{I}} - \frac {{\boldsymbol \nu } {\textbf{E}}}{\bar {\nu }}\biggr ) {\boldsymbol \nu } {\boldsymbol \eta }, \end{equation}

we can rewrite Eq. (4.10) in a more compact form as

(4.12) \begin{equation} {\mathrm{det}}(\sigma {\textbf{M}}_1 - {\textbf{B}}_1) = 0. \end{equation}

In the small-target limit, one has $\nu _j \ll 1$ so that ${\textbf{M}}_1$ is invertible since it is a small perturbation of the identity matrix. As a consequence, we get that

(4.13) \begin{equation} {\mathrm{det}}(\sigma {\textbf{I}} - {\textbf{M}}_1^{-1} {\textbf{B}}_1) = 0. \end{equation}

We now prove that the eigenvalues $\hat {\sigma }_j$ for $j=0,\ldots ,N-1$ of the matrix ${\textbf{M}}_1^{-1} {\textbf{B}}_1$ are real. To do so, we first write ${\textbf{M}}_1$ in Eq. (4.11) as

(4.14) \begin{equation} {\textbf{M}}_1 = {\textbf{I}} - \hat {{\textbf{B}}}_1 \hat {{\textbf{G}}} , \,\, \end{equation}

where we define the symmetric matrices $\hat {{\textbf{B}}}_1$ and $\hat {{\textbf{G}}}$ by

(4.15) \begin{equation} \hat {{\textbf{B}}}_1={\textbf{B}}_1{\boldsymbol \eta }^{-1} = \left ({\textbf{I}} - \frac {{\boldsymbol \nu } {\textbf{E}}}{\bar {\nu }}\right ){\boldsymbol \nu } , \qquad \hat {{\textbf{G}}}=-{\textbf{G}} - C_1 {\textbf{I}} . \end{equation}

By using the Neumann series to calculate ${\textbf{M}}_1^{-1}$ , which converges since $\nu _j\ll 1$ , we obtain that

(4.16) \begin{equation} {\textbf{M}}_1^{-1} \hat {{\textbf{B}}}_1 = \hat {{\textbf{B}}}_1 + \hat {{\textbf{B}}}_1\left (\sum _{n=1}^{\infty } {\textbf{K}}_n\right ) \hat {{\textbf{B}}}_1 , \qquad {\textbf{K}}_n = \left ( \hat {{\textbf{G}}} \hat {{\textbf{B}}}_1\right )^{n-1} \hat {{\textbf{G}}} . \end{equation}

Since ${\textbf{K}}_n$ is symmetric for each $n=1,2,\ldots$ , it follows that ${\textbf{M}}_1^{-1}\hat {{\textbf{B}}}_1$ is symmetric. Finally, we denote ${\textbf{D}}={\textbf{M}}_1^{-1}{\textbf{B}}_1 ={\textbf{M}}_{1}^{-1}\hat {{\textbf{B}}}_1 {\boldsymbol \eta }$ and introduce $\hat {{\textbf{D}}} = {\boldsymbol \eta }^{\frac 12} {\textbf{D}} {\boldsymbol \eta }^{-\frac 12} = {\boldsymbol \eta }^{\frac 12} {\textbf{M}}_1^{-1} \hat {{\textbf{B}}}_1 {\boldsymbol \eta }^{\frac 12}$ , which is symmetric, so that its eigenvalues are real. Since $\boldsymbol \eta$ is positive definite, $\hat {{\textbf{D}}}$ is related to $\textbf{D}$ by a similarity transformation so that the latter must also have real eigenvalues $\hat {\sigma }_j$ for $j=0,\ldots ,N-1$ .

Moreover, we note that $\hat {\sigma }_0 = 0$ . This follows since ${\textbf{B}}_1^\dagger {\textbf{e}} = 0$ , so that $\textbf{e}$ is an eigenvector of ${\textbf{B}}_1^\dagger$ and thus of $[{\textbf{M}}_1^\dagger ]^{-1} {\textbf{B}}_1^\dagger$ . The associated eigenvalue $0$ is thus an eigenvalue of ${\textbf{B}}_1 {\textbf{M}}_1^{-1}$ as well as of ${\textbf{M}}_1^{-1} {\textbf{B}}_1$ .

In summary, the eigenvalues of the matrix ${\textbf{M}}_1^{-1} {\textbf{B}}_1$ provide the leading terms in the asymptotic behaviour of the first $N$ eigenvalues $\sigma _j$ of the Steklov-Neumann problem:

(4.17) \begin{equation} \sigma _j \approx \hat {\sigma }_j, \qquad j\in \lbrace {1,\ldots ,N-1\rbrace }. \end{equation}

This is the main result of this section. We observe that $\sigma _0 = \hat {\sigma }_0 = 0$ , as expected.

4.3. Associated eigenfunctions

In addition, we construct the associated eigenfunctions of the first $N$ eigenvalues. For this purpose, let us rewrite Eq. (4.8) as

(4.18) \begin{equation} (\sigma {\textbf{I}} - {\textbf{M}}_1^{-1} {\textbf{B}}_1) {\textbf{A}} = \boldsymbol{0}. \end{equation}

While each eigenvalue $\hat {\sigma }_j$ of the matrix ${\textbf{M}}_1^{-1} {\textbf{B}}_1$ approximates the $j$ -th eigenvalue of the mixed Steklov-Neumann problem, the corresponding eigenvector of this matrix is the vector of coefficients $A_i$ determining the associated eigenfunction $V_j$ via Eq. (4.4), up to a multiplicative factor; we recall that $\chi$ is given by Eq. (4.7).

The missing multiplicative factor can be fixed by imposing an appropriate normalisation of eigenfunctions. For the Steklov problem, the natural normalisation is

(4.19) \begin{equation} \int \limits _{\Gamma } V^2\, ds = 1 , \end{equation}

where $\Gamma = \Gamma _{\varepsilon _1}\cup \cdots \cup \Gamma _{\varepsilon _N}$ . For the trivial eigenvalue $\sigma _0 = 0$ , one has a constant eigenfunction $V_0$ , whose normalisation yields: $V_0^2 = 1/|\Gamma |$ . In the following, we assume that $\sigma \gt 0$ .

To proceed, we recall that the inner solution near the $i$ -th patch reads in local coordinates is

(4.20) \begin{equation} V({{\boldsymbol{x}}}_i + \varepsilon _i \textbf{Q}_i^\dagger {\boldsymbol{y}}) \approx A_i g_{-\sigma \varepsilon _i}(\,{\boldsymbol{y}}). \end{equation}

Using the representation (3.4) of the Green’s function $g_{\mu }(\,{\boldsymbol{y}})$ , the restriction of $V$ onto $\Gamma _i$ becomes

(4.21) \begin{equation} V|_{\Gamma _i}(y_1) = V({{\boldsymbol{x}}}_i + \varepsilon _i \textbf{Q}_i^\dagger (y_1,0)^\dagger ) \approx \pi A_i \sum \limits _{k=0}^\infty \frac {\Psi _{2k}(\infty ) \Psi _{2k}(y_1,0)}{\mu _{2k} - \sigma \varepsilon _i} . \end{equation}

This equation helps to deduce the required condition on the coefficients $A_i$ :

\begin{align*} 1 & = \sum \limits _{i=1}^N \int \limits _{\Gamma _i} V^2 \, ds = \sum \limits _{i=1}^N \varepsilon _i \pi ^2 A_i^2 \sum \limits _{k=0}^\infty \frac {[\Psi _{2k}(\infty )]^2}{(\mu _{2k} - \sigma \varepsilon _i)^2} , \end{align*}

where we used the orthogonality of the eigenfunctions $\Psi _k$ (see Appendix C). The last sum can be re-written as the derivative of ${\mathcal{C}}(\mu )$ , denoted as ${\mathcal{C}}^{\prime }(\mu )$ :

(4.22) \begin{equation} 1 = - \pi \sum \limits _{i=1}^N \varepsilon _i A_i^2 {\mathcal{C}}^{\prime }({-}\sigma \varepsilon _i) . \end{equation}

When $0 \lt \sigma \varepsilon _i \ll 1$ , the Taylor expansion (3.8) implies ${\mathcal{C}}^{\prime }(\mu ) \approx -\pi /(2\mu ^2)$ and thus

(4.23) \begin{equation} \frac {2\sigma ^2}{\pi ^2} \approx \sum \limits _{i=1}^N \frac {A_i^2}{\varepsilon _i } . \end{equation}

To complete this section, let us briefly discuss the positivity of Steklov eigenfunctions on patches $\Gamma _{\varepsilon _j}$ . For a single patch, all Steklov eigenfunctions $V_j$ must change sign on the patch due to their orthogonality to $V_0 = 1/\sqrt {|\Gamma |}$ . When there are $N$ Steklov patches, the orthogonality still holds so that any eigenfunction $V_j$ with $j \gt 0$ must change sign on the union of patches $\Gamma = \Gamma _{\varepsilon _1} \cup \cdots \cup \Gamma _{\varepsilon _N}$ . However, it is generally unknown whether $V_j$ changes the sign or not on each patch $\Gamma _i$ . Looking at Eq. (4.21), one can expect that if $\sigma \varepsilon _i$ is small enough, the eigenfunction $V$ does not change sign on the patch $\Gamma _i$ (i.e., it is either positive or negative on it). Indeed, the term $1/({-}2\sigma \varepsilon _i)$ of the sum in Eq. (4.21) that corresponds to $k = 0$ , is expected to provide the dominant contribution as compared to the remaining terms. This property follows from the conjectured inequality (C.9). In other words, if the patches are small enough, the first $N$ eigenfunctions do not change their signs on each patch. This conjecture is confirmed by several numerical examples (not shown).

4.4. Example of two patches

When $N = 2$ , we calculate that

(4.24) \begin{equation} ({\textbf{I}} - {\boldsymbol \nu } {\textbf{E}}/\bar {\nu }){\boldsymbol \nu } = \gamma \left (\begin{array}{c@{\quad}c} 1 & -1 \\ -1 & 1 \\ \end{array} \right ), \end{equation}

where we label $\gamma = \nu _1\nu _2/(\nu _1+ \nu _2)$ . Then, from Eq. (4.11), we get

\begin{equation*} {\textbf{M}}_1 = {\textbf{I}} + \gamma \left (\begin{array}{c@{\quad}c} (R_1 - G_{1,2}) + C_1 & (G_{1,2} - R_2)-C_1 \\[4pt] (G_{1,2}-R_1) - C_1 & (R_2 - G_{1,2}) + C_1 \\ \end{array} \right ) , \qquad {\textbf{B}}_1 = \frac {\pi }{2} \gamma \left (\begin{array}{c@{\quad}c} 1/\varepsilon _1 & -1/\varepsilon _2 \\[4pt] -1/\varepsilon _1 & 1/\varepsilon _2 \\ \end{array} \right ). \end{equation*}

Since $\textbf{e}$ is a left-eigenvector of ${\textbf{M}}_1$ with eigenvalue one, the second eigenvalue of ${\textbf{M}}_1$ is simply $\mathrm{trace}({\textbf{M}}_1)-1$ . As a consequence, we find

(4.25) \begin{equation} {\mathrm{det}}({\textbf{M}}_1) = \mathrm{trace}({\textbf{M}}_1)-1=1 + \gamma [(R_1 + R_2 - 2G_{1,2}) + 2C_1] , \end{equation}

and

(4.26) \begin{equation} {\textbf{M}}_1^{-1} = \frac {1}{{\mathrm{det}}({\textbf{M}}_1)}\biggl [{\textbf{I}} + \gamma C_1 {\textbf{E}} + \gamma \left (\begin{array}{c@{\quad}c} R_2 - G_{1,2} & R_2-G_{1,2} \\[4pt] R_1 - G_{1,2} & R_1 - G_{1,2} \\ \end{array} \right )\biggr ] , \end{equation}

from which we calculate

(4.27) \begin{equation} {\textbf{M}}_1^{-1} {\textbf{B}}_1 = \frac {\pi \gamma /2}{{\mathrm{det}}({\textbf{M}}_1)} \left (\begin{array}{cc} 1/\varepsilon _1 & - 1/\varepsilon _2 \\[4pt] -1/\varepsilon _1 & 1/\varepsilon _2 \\ \end{array} \right ). \end{equation}

The two eigenvalues of this matrix are $\hat {\sigma }_0 = 0$ and

(4.28) \begin{equation} \hat {\sigma }_1 = \frac {\pi \gamma (1/\varepsilon _1 + 1/\varepsilon _2)}{2\,{\mathrm{det}}({\textbf{M}}_1)} , \end{equation}

so that upon solving for $\hat {\sigma }_1^{-1}$ , we get

(4.29) \begin{equation} \frac {1}{\hat {\sigma }_1} = \frac {2\varepsilon _1\varepsilon _2}{\pi (\varepsilon _1+\varepsilon _2)} \biggl [\frac {1}{\gamma } + (R_1 + R_2 - 2G_{1,2}) + 2C_1\biggr ] . \end{equation}

To simplify this expression, we use

(4.30) \begin{equation} \gamma = \frac {1}{1/\nu _1 + 1/\nu _2} = \frac {1}{-\ln (\varepsilon _1\varepsilon _2)} . \end{equation}

This yields the following asymptotic behaviour for the Steklov eigenvalue:

(4.31) \begin{align} \frac {1}{\sigma _1} \approx \frac {2\varepsilon _1\varepsilon _2}{\pi (\varepsilon _1+\varepsilon _2)} [{-}\ln (\varepsilon _1\varepsilon _2) + 2C_1 + (R_1 + R_2 - 2G_{1,2})]. \end{align}

For instance, if $\Omega$ is the unit disk, Eq. (2.28) yields

(4.32) \begin{align} \frac {1}{\sigma _1} \approx \frac {2\varepsilon _1\varepsilon _2}{\pi (\varepsilon _1+\varepsilon _2)} ({-}\ln (\varepsilon _1\varepsilon _2) + 2C_1 + 2\ln |{{\boldsymbol{x}}}_1-{{\boldsymbol{x}}}_2|). \end{align}

Moreover, the corresponding eigenfunction $V_1$ can be easily found by noting that $A_1 = -A_2$ from Eq. (2.11), whereas Eq. (4.23) implies

(4.33) \begin{equation} A_1 = -A_2 \approx \frac {\sqrt {2} \sigma _1}{\pi \sqrt {1/\varepsilon _1 + 1/\varepsilon _2}} . \end{equation}

These coefficients determine $(V_1)|_{\Gamma _i}$ via Eq. (4.21).

4.5. Example of identical equally spaced patches on the boundary of the unit disk

When all patches are of the same size, $\varepsilon _j = \varepsilon$ , we have ${\boldsymbol \nu } = \nu {\textbf{I}}$ and ${\boldsymbol \eta } = {\textbf{I}} \pi /(2\varepsilon )$ , so that

(4.34) \begin{equation} {\textbf{M}}_1 = {\textbf{I}} + \nu \biggl ({\textbf{I}} - \frac {{\textbf{E}}}{N}\biggr ) ({\textbf{G}} + C_1 {\textbf{I}}), \qquad {\textbf{B}}_1 = \frac {\pi \nu }{2\varepsilon } \biggl ({\textbf{I}} - \frac {{\textbf{E}}}{N}\biggr ) . \end{equation}

If the patches are equally spaced on the boundary of the unit disk, $\textbf{G}$ is circulant and symmetric, and its eigenvectors and eigenvalues were given in Eqs. (2.32, 2.33). As shown earlier in Section 2.5, one has $({\textbf{I}} - {\textbf{E}}/N) \textbf{q}_j = \textbf{q}_j$ for any $j = 1,2,\ldots ,N-1$ , and $({\textbf{I}} - {\textbf{E}}/N) \textbf{q}_N = 0$ . As a consequence, ${\textbf{M}}_1 \textbf{q}_N = \textbf{q}_N$ and ${\textbf{M}}_1 \textbf{q}_j = (1 + \nu (\kappa _j+C_1)) \textbf{q}_j$ , where $\kappa _j$ are given explicitly by Eq. (2.40). Upon calculating ${\textbf{M}}_1^{-1}\textbf{q}_j$ , we readily find that

\begin{equation*} {\textbf{B}}_1 {\textbf{M}}_1^{-1} \textbf{q}_N = 0, \qquad {\textbf{B}}_1 {\textbf{M}}_1^{-1} \textbf{q}_j = \frac {\pi \nu }{2\varepsilon } \bigl (1 + \nu (\kappa _j+C_1)\bigr )^{-1} \textbf{q}_j . \end{equation*}

Since the eigenvalues of the matrices ${\textbf{B}}_1 {\textbf{M}}_1^{-1}$ and ${\textbf{M}}_1^{-1}{\textbf{B}}_1$ are identical, we have from Eq. (4.13) that $\hat {\sigma }_0 = 0$ and

(4.35) \begin{equation} \hat {\sigma }_j = \frac {\pi \nu }{2\varepsilon [1 + \nu (\kappa _j+C_1)]}, \qquad j\in \lbrace {1,\ldots ,N-1\rbrace }. \end{equation}

Substituting $C_1={3/2}-\ln {2}$ from Eq. (3.9) and $\nu ={-1/\ln \varepsilon }$ into Eq. (4.35), we obtain the following small- $\varepsilon$ asymptotic result for the first $N-1$ eigenvalues of the mixed Steklov-Neumann problem:

(4.36) \begin{equation} \frac {1}{\varepsilon \sigma _j} \approx \frac {2}{\pi } \biggl ({-}\ln (\varepsilon ) + \frac 32 - \ln (2) + \kappa _j \biggr ), \end{equation}

for $j = 1,2,\ldots ,N-1$ . When the number of patches is large, i.e. $N \gg 1$ , while still enforcing the well-separated patch assumption $\varepsilon N\ll \pi$ , we can use the asymptotic relation (B.12) for $\kappa _j$ , valid for $j\ll N$ , to conclude that

(4.37) \begin{equation} \frac {1}{\varepsilon \sigma _j} \approx \frac {2}{\pi } \biggl ( \frac {N}{2j} -\ln \left (\frac {N\varepsilon }{\pi }\right ) - \frac {1}{3}\biggr ). \end{equation}

If $N\varepsilon$ is not too small, the logarithmic and constant terms can be neglected to yield to a first approximation

(4.38) \begin{equation} \sigma _j \approx \frac {\pi j}{N\varepsilon } \quad \mbox{for} \quad j\ll N, \end{equation}

when $N\gg 1$ . It is instructive to compare this approximation to the case of a single Steklov patch of half-length $\varepsilon _1 = N\varepsilon$ , for which $\sigma _j \approx \mu _j/\varepsilon _1 \approx (\pi j/2)/(N\varepsilon )$ , where we used the asymptotic relation (C.11). As a consequence, the configuration with a single patch of half-length $N\varepsilon$ yields approximately twice smaller eigenvalues. This suggests that the fragmentation of a patch will increase the eigenvalues.

5.1. Matched asymptotic analysis

Expectedly, we can combine formerly derived inner solutions for Dirichlet and Steklov patches, whereas the outer solution is still written as the linear combination (4.4). As a consequence, we must enforce the singularity behaviour

(5.3) \begin{equation} V \sim A_j \left[\ln |{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_j| + 1/\nu _j + \delta _{j,1} {\mathcal{C}}({-}\sigma \varepsilon _1) + o(1)\right] , \end{equation}

as ${{\boldsymbol{x}}}\to {{\boldsymbol{x}}}_j$ for all patches $j\in \lbrace {1,\ldots ,N\rbrace }$ . Here, we have used the former definition $\nu _j = -1/\ln (\varepsilon _j/2)$ for Dirichlet patches and $\nu _1 = -1/\ln (\varepsilon _1)$ for the Steklov patch. By ensuring that $V$ in Eq. (4.4) satisfies Eq. (5.3), we obtain that

(5.4) \begin{equation} A_j + \nu _j \left[A_j R_j + \sum \limits _{i=1 \atop i\ne j}^N G_{j,i} A_i\right] = \chi \nu _j - \nu _1 A_1 {\mathcal{C}}({-}\sigma \varepsilon _1) \delta _{j,1} . \end{equation}

Together with Eq. (2.11), they form a system of $N+1$ linear equations for the unknowns $A_i$ and $\chi$ . Using the former matrix notation in Eq. (2.15), Eq. (5.4) becomes

(5.5) \begin{equation} \bigl ({\textbf{I}} + {\boldsymbol \nu } {\textbf{G}}\bigr ) {\textbf{A}} = \chi {\boldsymbol \nu } {\textbf{e}} - \nu _1 {\mathcal{C}}({-}\sigma \varepsilon _1) {\textbf{e}}_1 {\textbf{e}}_1^\dagger {\textbf{A}}, \end{equation}

with ${\textbf{e}}^\dagger {\textbf{A}} = 0$ . Upon left-multiplying by ${\textbf{e}}^\dagger$ , we get

(5.6) \begin{equation} \chi = \frac {1}{\bar {\nu }} [{\textbf{e}}^\dagger {\boldsymbol \nu } {\textbf{G}} {\textbf{A}} + {\mathcal{C}}({-}\sigma \varepsilon _1) \nu _1 {\textbf{e}}^\dagger {\textbf{E}}_1 {\textbf{A}}], \end{equation}

where $\bar {\nu }$ was defined in Eq. (2.18), and we introduced the matrix ${\textbf{E}}_1 = {\textbf{e}}_1 {\textbf{e}}_1^\dagger$ for a shorter notation. Eliminating $\chi$ from Eq. (5.5), we obtain that

(5.7) \begin{equation} \biggl ({\textbf{I}} + {\boldsymbol \nu } \biggl ({\textbf{I}} - \frac {{\textbf{E}} {\boldsymbol \nu }}{\bar {\nu }} \biggr ) {\textbf{G}} \biggr ) {\textbf{A}} + {\mathcal{C}}({-}\sigma \varepsilon _1) \nu _1 \biggl ({\textbf{I}} - \frac {{\boldsymbol \nu } {\textbf{E}}}{\bar {\nu }} \biggr ) {\textbf{E}}_1 {\textbf{A}} = \boldsymbol{0} , \end{equation}

where we recall that ${\textbf{E}} = {\textbf{e}} {\textbf{e}}^\dagger$ . Upon introducing $\textbf{B}$ by

(5.8) \begin{equation} {\textbf{B}} = \biggl ({\textbf{I}} - \frac {{\boldsymbol \nu } {\textbf{E}}}{\bar {\nu }} \biggr ) {\textbf{E}}_1 , \end{equation}

and using the matrix ${\textbf{M}}_0$ from Eq. (2.20), we rewrite the matrix system in Eq. (5.7) as

(5.9) \begin{equation} \bigl ({\textbf{M}}_0 + {\mathcal{C}}({-}\sigma \varepsilon _1) \nu _1 {\textbf{B}}\bigr ) {\textbf{A}} = \boldsymbol{0} . \end{equation}

The condition, under which this matrix equation admits a nontrivial solution is

(5.10) \begin{equation} {\mathrm{det}}\bigl ({\textbf{M}}_0 + {\mathcal{C}}({-}\sigma \varepsilon _1) \nu _1 {\textbf{B}}\bigr ) = 0, \end{equation}

which is a scalar problem that determines $\sigma$ .

To rewrite this problem in a more explicit form, we first observe that ${\mathrm{rank}}({\textbf{B}}) = 1$ since ${\textbf{E}}_1 \textbf{q} = 0$ for any vector $\textbf{q} \in {\mathbb{R}}^N$ such that ${\textbf{e}}_1^\dagger \textbf{q} = 0$ . Since ${\textbf{e}}^\dagger {\textbf{e}}_1 = 1$ , we can rewrite $\textbf{B}$ in a more convenient rank-one form as

(5.11) \begin{equation} {\textbf{B}} = \biggl ({\textbf{e}}_1 - \frac {{\boldsymbol \nu } {\textbf{e}}}{\bar {\nu }} \biggr ) {\textbf{e}}_1^\dagger = {\textbf{a}} {\textbf{b}}^\dagger , \quad \mbox{where} \quad {\textbf{a}} = {\textbf{e}}_1 - \frac {{\boldsymbol \nu } {\textbf{e}}}{\bar {\nu }} , \qquad {\textbf{b}} = {\mathcal{C}}({-}\sigma \varepsilon _1) \nu _1 {\textbf{e}}_1 . \end{equation}

To proceed, we need the matrix determinant lemma [Reference Ding and Zhu23].

Lemma. Let ${\textbf{M}} = {\textbf{M}}_0 + {\textbf{a}} {\textbf{b}}^\dagger$ be a perturbation of an invertible matrix ${\textbf{M}}_0 \in {\mathbb{R}}^{N,N}$ by a rank-one matrix ${\textbf{a}} {\textbf{b}}^\dagger$ . Then

(5.12) \begin{equation} {\mathrm{det}}({\textbf{M}}_0 + {\textbf{a}} {\textbf{b}}^\dagger ) = (1 + {\textbf{b}}^\dagger {\textbf{M}}_0^{-1} {\textbf{a}})\, {\mathrm{det}}({\textbf{M}}_0). \end{equation}

It follows that ${\mathrm{det}}({\textbf{M}}) = 0$ if and only if ${\textbf{b}}^\dagger {\textbf{M}}_0^{-1} {\textbf{a}} = -1$ .

In the small-target limit, all $\nu _j \ll 1$ so that the matrix ${\textbf{M}}_0$ is invertible since it is a small perturbation of the identity matrix in Eq. (2.20). Applying the lemma above to our setting, we determine the condition on $\sigma$ as

(5.13) \begin{equation} {\mathcal{C}}({-}\sigma \varepsilon _1) = C , \end{equation}

where

(5.14) \begin{equation} C = - \frac {1}{\nu _1} \biggl ( {\textbf{e}}_1^\dagger {\textbf{M}}_0^{-1} \biggl [{\textbf{e}}_1 - \frac {{\boldsymbol \nu } {\textbf{e}}}{\bar {\nu }}\biggr ]\biggr )^{-1} , \end{equation}

with the vectors and matrices $\textbf{e}$ , ${\textbf{e}}_1$ , $\boldsymbol \nu$ , and ${\textbf{M}}_0$ being defined in Eqs. (2.15, 2.20). This is the main result of this section that will allow us to determine the asymptotic behaviour of the Steklov eigenvalues and their dependence on the configuration and sizes of all patches that are captured via the constant $C$ in Eq. (5.14). We further emphasise that the homogeneous matrix equation (5.9) cannot uniquely determine the coefficients $A_i$ . In fact, an eigenfunction $V$ can be found up to a multiplicative factor that has to be fixed by normalisation (see below).

As stated above, the matrix ${\textbf{M}}_0$ is a small perturbation of the identity matrix in the small-target limit, so that ${\textbf{M}}_0^{-1} \sim {\textbf{I}}$ to leading order, which implies that ${\textbf{e}}_1^\dagger {\textbf{M}}_0^{-1} [{\textbf{e}}_1 - {\boldsymbol \nu } {\textbf{e}}/\bar {\nu }] \sim 1 - \nu _1/\bar {\nu } \gt 0$ . Since $0 \lt \nu _1 \ll 1$ , we conclude that the constant $C$ in Eq. (5.14) is negative and large:

(5.15) \begin{equation} C \lt 0, \qquad |C| \gg 1. \end{equation}

5.2. Asymptotic behaviour of eigenvalues and eigenfunctions

Denoting $\mu = -\sigma \varepsilon _1$ , we recast Eq. (5.13) as

(5.16) \begin{equation} {\mathcal{C}}(\mu ) = C . \end{equation}

The spectral expansion (3.7) of the function ${\mathcal{C}}(\mu )$ allows one to solve this equation numerically for any fixed negative value $C$ given by Eq. (5.14). Since the derivative ${\mathcal{C}}^{\prime }(\mu ) = d{\mathcal{C}}(\mu )/d\mu$ is negative, ${\mathcal{C}}(\mu )$ is a continuous and monotonically decreasing function on each interval $({-}\mu _{2j+2}, -\mu _{2j})$ , with $j = 0,1,\ldots$ . Moreover, it ranges from $+\infty$ to $-\infty$ on each interval. As a consequence, for any fixed value $C$ , there exist infinitely many negative solutions of Eq. (5.16), denoted as $-\hat {\mu }_{2j}$ , such that

(5.17) \begin{equation} \mu _{2j} \leq \hat {\mu }_{2j} \leq \mu _{2j+2} \quad \mbox{for} \, j = 0, 1, \ldots . \end{equation}

This property facilitates the numerical solution, as a single zero has to be searched on each interval. Moreover, as the coefficients $[\Psi _{2k}(\infty )]^2$ are small (see Table C1) and decrease with $k$ , whereas $C$ is negative and large, one has $\hat {\mu }_{2j} \approx \mu _{2j}$ for $j \gt 0$ . This is consistent with the intuitive picture that large Steklov eigenvalues become insensitive to Dirichlet patches in the small-target limit, and one retrieves the asymptotic behaviour for a single Steklov patch [Reference Grebenkov38].

The solutions $\hat {\mu }_{2j} \approx \mu _{2j}$ determine the leading-order term in the asymptotic behaviour of the eigenvalues $\sigma _{2j}$ :

(5.18) \begin{equation} \sigma _{2j} \approx \frac {\mu _{2j}}{\varepsilon _1} \qquad \mbox{for} \, j = 1,2,\ldots . \end{equation}

In contrast, the smallest eigenvalue $\sigma _0$ involves the solution $\hat {\mu }_0$ , which may actually be small in the small-target limit. We discuss this case separately in Section 5.3.

We also mention that the analysis above provides the leading-order approximation to the associated Steklov eigenfunction, restricted to $\Gamma _{\varepsilon _1}$ . We recall that the inner solution near the Steklov patch is $V_{2j} \sim A_1 g_{-\sigma _{2j} \varepsilon _1}(\,{\boldsymbol{y}})$ , with the Green’s function $g_\mu (\,{\boldsymbol{y}})$ given by Eq. (3.4). As a consequence, its restriction onto the patch reads

(5.19) \begin{equation} V_{2j}\biggr |_{\Gamma _{\varepsilon _1}} \approx a_{2j} \sum \limits _{k=0}^\infty \frac {\Psi _{2k}(\infty )}{\mu _{2k} - \varepsilon _1 \sigma _{2j}} \Psi _{2k}(\,{\boldsymbol{y}}), \end{equation}

where the proportionality coefficient $a_{2j}$ is fixed by the conventional normalisation of the Steklov eigenfunction:

(5.20) \begin{equation} 1 = \int \limits _{\Gamma _{\varepsilon _1}} V_{2j}^2\, ds \approx \varepsilon _1 a_{2j}^2 \sum \limits _{k=0}^\infty \frac {[\Psi _{2k}(\infty )]^2}{(\mu _{2k} - \varepsilon _1 \sigma _{2j})^2} , \end{equation}

where the orthogonality of $\Psi _{2k}$ was used. Since $\varepsilon _1 \sigma _{2j} \approx \mu _{2j}$ , the eigenfunction $\Psi _{2j}$ provides the dominant contribution, and one gets for each $j\in \lbrace {1,2,\ldots \rbrace }$ that

(5.21) \begin{equation} V_{2j}({{\boldsymbol{x}}}_1 + \varepsilon _1 \textbf{Q}_1^\dagger (y_1,0)^\dagger ) \approx \frac {1}{\sqrt {\varepsilon _1}} \Psi _{2j}(y_1,0) \end{equation}

on the Steklov patch (i.e., for $|y_1| \leq 1$ ).

We emphasise that the analysis above allowed us to access only half of eigenvalues with even indices $2j$ that correspond to symmetric eigenmodes. In turn, the eigenvalues with odd indices correspond to antisymmetric eigenmodes, for which the integral over the Steklov patch is zero. As discussed at the beginning of Section 5, such eigenfunctions vanish away from the Steklov patch so that their asymptotic behaviour can be determined directly from the local solution. As a consequence, we get a leading-order approximation

(5.22) \begin{equation} \sigma _{2j+1} \approx \frac {\mu _{2j+1}}{\varepsilon _1} \qquad \mbox{for} \, j=0,1,\ldots , \end{equation}

and on the Steklov patch we have

(5.23) \begin{equation} V_{2j+1}({{\boldsymbol{x}}}_1 + \varepsilon _1 \textbf{Q}_1^\dagger (y_1,0)^\dagger ) \approx \frac {1}{\sqrt {\varepsilon _1}} \Psi _{2j+1}(y_1,0) . \end{equation}

In summary, our analysis justifies theoretically the intuitively appealing scaling argument that the eigenvalues $\sigma _j$ and eigenfunctions $(V_j)|_{\Gamma _{\varepsilon _1}}$ (on the Steklov patch) can be approximated by the eigenvalues $\mu _j$ and eigenfunctions $\Psi _j$ of the auxiliary problem (3.5), as if there were no Dirichlet patches. Regardless of the symmetry of eigenfunctions (and parity of its index), we can combine the former leading-order approximations as

(5.24a) \begin{align} \sigma _j & \approx \frac {\mu _j}{\varepsilon _1} \quad \mbox{for} \, j=1,2,\ldots , \\[-5pt] \nonumber \end{align}

(5.24b) \begin{align} V_j({{\boldsymbol{x}}}_1 + \varepsilon _1 \textbf{Q}_1^\dagger (y_1,0)^\dagger ) & \approx \frac {1}{\sqrt {\varepsilon _1}} \Psi _j(y_1,0) . \\[8pt] \nonumber \end{align}

In contrast, the presence of Dirichlet patches must affect the principal eigenvalue $\sigma _0$ and the associated eigenfunction $V_0$ , as explained below.

5.3. The principal eigenvalue

Since $\mu _0 = 0$ , the smallest solution of Eq. (5.16), $\hat {\mu }_0$ , is close to $0$ . Indeed, as the constant $C$ is large and negative, one needs to have $|\mu | \ll 1$ to ensure that ${\mathcal{C}}(\mu )$ is also large and negative. Under the condition $|\mu |\ll 1$ , we can use the approximation (3.10), which can be easily inverted to get the explicit result

(5.25) \begin{equation} \frac {1}{\mu } \approx \frac {2}{\pi } \bigl [{\mathcal{C}}(\mu ) - C_1 \bigr ] = \frac {2}{\pi } \bigl [C - C_1\bigr ] . \end{equation}

As a consequence, substitution of Eq. (5.14) here yields the asymptotic behaviour of the principal eigenvalue $\sigma _0$ :

(5.26) \begin{equation} \frac {1}{\varepsilon _1 \sigma _0} \approx - \frac {2}{\pi } \biggl [\ln (\varepsilon _1) \biggl ( {\textbf{e}}_1^\dagger {\textbf{M}}_0^{-1} \biggl [{\textbf{e}}_1 - \frac {{\boldsymbol \nu } {\textbf{e}}}{\bar {\nu }}\biggr ]\biggr )^{-1} - C_1\biggr ]. \end{equation}

In sharp contrast to Eq. (5.18) for $\sigma _j$ with $j \geq 1$ , the principal eigenvalue $\sigma _0$ exhibits a slower divergence ${\mathcal O}(1/(\varepsilon _1\ln (\varepsilon _1)))$ . This is one of the main results of this section. The associated eigenfunction $V_0$ is given by Eq. (5.19) with the normalisation condition (5.20). We stress that Eq. (5.19) cannot be reduced to the approximation (5.21) in this case.

5.4. Example of two patches

When there are two patches ( $N = 2$ ), Eqs. (5.14, 5.26) can be readily solved. Substituting ${\textbf{M}}_0^{-1}$ from Eq. (2.24) and

(5.27) \begin{equation} {\textbf{e}}_1 - \frac {{\boldsymbol \nu } {\textbf{e}}}{\bar {\nu }} = \frac {\nu _2}{\nu _1+\nu _2} \left (1,-1\right )^{\dagger } \end{equation}

into Eq. (5.14), we get after simplifications that

(5.28) \begin{equation} C = \ln (\varepsilon _1 \varepsilon _2/2) - (R_1+R_2-2G_{1,2}) . \end{equation}

Using Eq. (5.25) with $C_1={3/2}-\ln {2}$ from Eq. (3.9), the asymptotic behaviour of the principal eigenvalue is

(5.29) \begin{equation} \frac {1}{\varepsilon _1 \sigma _0} \approx \frac {2}{\pi }\biggl ( -\ln (\varepsilon _1 \varepsilon _2) + \frac 32 + (R_1+R_2-2G_{1,2})\biggr ) . \end{equation}

For instance, when $\Omega$ is the unit disk, one can substitute Eq. (2.28) into Eq. (5.29) to get

(5.30) \begin{equation} \frac {1}{\varepsilon _1 \sigma _0} \approx \frac {2}{\pi } \biggl ( -\ln (\varepsilon _1 \varepsilon _2) + \frac 32 + 2\ln |{{\boldsymbol{x}}}_1 - {{\boldsymbol{x}}}_2| \biggr ) . \end{equation}

Figure 7 illustrates the remarkable accuracy of this asymptotic relation.

Figure 8 shows the behaviour of the Steklov eigenfunctions $V_j$ restricted onto the Steklov patch. For the principal eigenmode with $j = 0$ , this restriction is positive, as expected. The asymptotic formula (5.19) yields an accurate approximation. Let us now look at other eigenmodes with $j = 1,2,3$ , for which $\sigma _j \approx \mu _j/\varepsilon _1$ . We see that the restriction of $V_j$ and its approximation (5.24) are in excellent agreement, for both symmetric and antisymmetric eigenfunctions, even though both considered patches are not small.

5.5. Example of equally spaced patches on the unit disk

We now consider another setting where the general formulas (5.13) and (5.14) can be simplified. We suppose that all Dirichlet patches with $j = 2,\ldots ,N$ are of the same length, so that $\varepsilon _j = \varepsilon$ for $j\in \lbrace {2,\ldots ,N\rbrace }$ , whereas for the Steklov patch we have $\varepsilon _1 = \ell _1 \varepsilon /2$ , for some $\ell _1\gt 0$ . To treat this case, we can impose in our general formula (5.14) for $C$ that

(5.31) \begin{equation} \nu =\nu _1 = \ldots = \nu _N = -{1/\ln (\varepsilon /2)}, \qquad {\boldsymbol \nu }=\nu {\textbf{I}} , \end{equation}

provided that we shift ${\mathcal{C}}({-}\sigma \varepsilon _1)$ appropriately in the relation (5.13). To determine this shift, we write the singularity condition (5.3) as ${{\boldsymbol{x}}}\to {{\boldsymbol{x}}}_1$ associated with the Steklov patch $j=1$ as

\begin{align*} V & \sim A_1 [\ln |{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_1| - \ln \varepsilon _1 + {\mathcal{C}}({-}\sigma \varepsilon _1) + o(1)] \\[4pt] & \sim A_1 [\ln |{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_1| + 1/\nu + \tilde {{\mathcal{C}}}({-}\sigma \varepsilon _1) + o(1)], \end{align*}

where we have defined $\tilde {{\mathcal{C}}}({-}\sigma \varepsilon _1)$ by

(5.32) \begin{equation} \tilde {{\mathcal{C}}}({-}\sigma \varepsilon _1)={\mathcal{C}}({-}\sigma \varepsilon _1) - \ln \left (\frac {2\varepsilon _1}{\varepsilon }\right ) . \end{equation}

As a result, by repeating the steps of Section 5.1, we need only replace Eq. (5.13) by

(5.33) \begin{equation} {\mathcal{C}}({-}\sigma \varepsilon _1) - \ln \left ( \frac {2\varepsilon _1}{\varepsilon } \right ) = C , \end{equation}

where $C$ is given by Eq. (5.14) with ${\boldsymbol \nu } = \nu {\textbf{I}}$ , $\bar {\nu } = N\nu$ , and $\nu ={-1/\ln \left ({\varepsilon /2}\right )}$ , which yields

(5.34) \begin{equation} C = \frac {1}{\nu } \biggl (\frac {{\textbf{e}}_1^\dagger {\textbf{M}}_0^{-1} {\textbf{e}}}{N} - {\textbf{e}}_1^\dagger {\textbf{M}}_0^{-1} {\textbf{e}}_1\biggr )^{-1} . \end{equation}

Figure 6. Illustration for the unit disk with Steklov and Dirichlet patches. (a) One Steklov patch of length $2\varepsilon _1 = 0.4$ at ${{\boldsymbol{x}}}_1 = (1,0)$ (blue) and one Dirichlet patch of length $2\varepsilon _2 = 0.6$ (red), whose centre ${{\boldsymbol{x}}}_2$ is at angle $\theta = 2\pi /3$ . (b) One Steklov patch of length $2\varepsilon _1 = 2\varepsilon = 0.2$ at ${{\boldsymbol{x}}}_1 = (1,0)$ (blue) and three Dirichlet patches of length $2\varepsilon _j = 0.4$ (red), whose centres ${{\boldsymbol{x}}}_j$ are equally spaced on the boundary of the unit disk.

Figure 7. Dependence of $1/(\varepsilon _1\sigma _0)$ on $\varepsilon _2$ for the unit disk with a Steklov patch of length $2\varepsilon _1$ (located at ${{\boldsymbol{x}}}_1 = (1,0)$ ), and one Dirichlet patch of length $2\varepsilon _2$ , located at ${{\boldsymbol{x}}}_2$ . Symbols present the numerical solution by a FEM with the maximal mesh size $h_{\textrm {max}} = 0.005$ and lines show Eq. (5.30). (a) ${{\boldsymbol{x}}}_2 = (0,1)$ and three values of $\varepsilon _1$ : $\varepsilon _1 = \pi /6$ (circles), $\varepsilon _1 = \pi /12$ (squares), and $\varepsilon _1 = \pi /24$ (triangles). (b) $\varepsilon _1 = \pi /12$ and two locations of the Dirichlet patch: ${{\boldsymbol{x}}}_2 = ({-}1,0)$ (circles, $\theta = \pi$ ), and ${{\boldsymbol{x}}}_2 = (0,1)$ (squares, angle $\theta = \pi /2$ ).

If all ${{\boldsymbol{x}}}_j$ are equally spaced on the boundary of the unit disk (such as shown in Figure 6(b)), then the matrix $\textbf{G}$ is circulant and symmetric, so that its eigenvectors and eigenvalues are known exactly (see Section 2.5). Moreover, the matrix ${\textbf{M}}_0$ admits a spectral representation and thus can be inverted explicitly. Using Eq. (2.36), we get

(5.35) \begin{equation} {\textbf{e}}_1^\dagger {\textbf{M}}_0^{-1} {\textbf{e}} = 1 , \qquad {\textbf{e}}_1^\dagger {\textbf{M}}_0^{-1} {\textbf{e}}_1 = ({\textbf{e}}_1^\dagger \textbf{q}_N) (\textbf{q}_N^\dagger {\textbf{e}}_1) + \sum \limits _{j=1}^{N-1} \frac {({\textbf{e}}_1^\dagger \textbf{q}_j) \, (\textbf{q}_j^\dagger {\textbf{e}}_1)}{1 + \nu \kappa _j} , \end{equation}

where $\textbf{q}_j$ and $\kappa _j$ were defined by Eqs. (2.32, 2.40). However, since $({\textbf{e}}_1^\dagger \textbf{q}_j) = \omega ^j/\sqrt {N}$ and $(\textbf{q}_1^\dagger {\textbf{e}}_1) = \omega ^{-j}/\sqrt {N}$ , we get

(5.36) \begin{equation} {\textbf{e}}_1^\dagger {\textbf{M}}_0^{-1} {\textbf{e}}_1 = \frac {1}{N}\left[1 + \sum \limits _{j=1}^{N-1} (1 + \nu \kappa _j)^{-1} \right]. \end{equation}

Substituting this expression together with Eq. (5.35) into Eq. (5.34), we find

(5.37) \begin{equation} C = - \frac {N}{\nu } \left(\sum \limits _{j=1}^{N-1} \frac {1}{1 + \nu \kappa _j} \right)^{-1} . \end{equation}

As a result, Eqs. (5.25, 5.33) with $C_1={3/2}-\ln {2}$ imply

(5.38) \begin{align} \frac {1}{\varepsilon _1 \sigma _0} & \approx - \frac {2}{\pi } [{\mathcal{C}}({-}\sigma \varepsilon _1) - C_1 ] = - \frac {2}{\pi } \left[C + \ln \left (\frac {2\varepsilon _1}{\varepsilon }\right ) - C_1\right]\nonumber \\[5pt] & \approx \frac {2}{\pi } \left[\frac {N}{\sum \nolimits _{j=1}^{N-1} \bigl (\ln (2/\varepsilon ) + \kappa _j\bigr )^{-1}} + C_1 - \ln \left (\frac {2\varepsilon _1}{\varepsilon }\right )\right]. \end{align}

Figure 9 illustrates the behaviour of $1/(\varepsilon _1\sigma _0)$ as a function of $\ln (\varepsilon )$ for the unit disk with one Steklov patch of length $2\varepsilon$ , and several Dirichlet patches of length $4\varepsilon$ , which are equally spaced on the boundary of the unit disk. We observe an excellent agreement between the asymptotic formula (5.38) and numerical results.

Figure 8. The eigenfunctions $V_j$ restricted on the Steklov patch $\Gamma _{\varepsilon _1}$ , for the unit disk with a Steklov patch of length $2\varepsilon _1 = \pi /12 \approx 0.26$ (located at ${{\boldsymbol{x}}}_1 = (1,0)$ ), and one Dirichlet patch of length $2\varepsilon _2 = \pi /6 \approx 0.52$ , located at ${{\boldsymbol{x}}}_2 = ({-}1,0)$ . Filled circles present the numerical solution by a FEM with the maximal meshsize $h_{\textrm {max}} = 0.005$ , while solid lines show Eq. (5.19) for $j = 0$ and Eq. (5.24) for $j \gt 0$ . Four panels present the cases $j = 0,1,2,3$ .

Figure 9. Dependence of $1/(\varepsilon \sigma _0)$ on $\varepsilon$ for the unit disk with one Steklov patch of length $2\varepsilon$ (located at ${{\boldsymbol{x}}}_0 = (1,0)$ ), and $N-1$ Dirichlet patches of length $4\varepsilon$ , equally spaced on the boundary of the unit disk. Symbols present the numerical solution by a FEM with the maximal meshsize $h_{\textrm {max}} = 0.005$ and lines show Eq. (5.38).

Figure 10. Dependence of $1/(\varepsilon _1 \sigma )$ on $\varepsilon$ for the unit disk with the Steklov patch of length $2\varepsilon _1 = 0.2$ and $63$ Dirichlet patches (each of length $2\varepsilon$ ) that are equally spaced on the boundary of the unit disk. Filled circles correspond to $\kappa _j$ obtained via the discrete sum (5.38), the solid line indicates its large- $N$ approximation (5.43), and the dashed line is the low-order approximation (5.44).

Although the eigenvalues $\kappa _j$ are known explicitly via Eq. (2.40), it is instructive to inspect their asymptotic behaviour for large $N$ (see Appendix B). We aim at approximating the sum in the denominator of Eq. (5.38):

(5.39) \begin{equation} S = \frac {1}{N} \sum \limits _{j=1}^{N-1} \frac {1}{\ln (2/\varepsilon ) + \kappa _j} . \end{equation}

The degeneracy $\kappa _{N-j} = \kappa _j$ allows us to limit this sum to $j \leq {N/2}$ when $N$ is even. By using the asymptotic result in Eq. (B.13) for $\kappa _j$ when $N\gg 1$ , we find that

(5.40) \begin{equation} \ln \left ({2/\varepsilon }\right ) + \kappa _j \approx \frac {1}{2\xi } + \zeta + a \frac {\pi ^2\xi ^2}{9} , \end{equation}

where we have defined $\xi ={j/N}$ ,

(5.41) \begin{equation} \zeta = -\ln \left ( \frac {N\varepsilon }{b}\right )\!, \quad \mbox{and} \quad b = 4\pi e^{-11/6}\approx 2.009. \end{equation}

Here, the coefficient $a = 1.25$ was empirically introduced to improve the accuracy of the approximation of $\kappa _j$ (see Appendix B for details). In terms of $\zeta$ , and substituting ${j/N} = \xi$ , we view $S$ as a Riemannian approximation of the integral defined by

(5.42) \begin{equation} S \approx S(\zeta ) = 4 \int \limits _0^{1/2} \frac {\xi \, d\xi }{1 + 2\xi \left (\zeta + a \pi ^2\xi ^2/9\right )} . \end{equation}

As a result, Eq. (5.38) is approximated for $N\gg 1$ by

(5.43) \begin{align} \frac {1}{\varepsilon _1 \sigma _0} & \approx \frac {2}{\pi } \biggl [ \frac {1}{S(\zeta )} + C_1 - \ln \left (\frac {2\varepsilon _1}{\varepsilon }\right )\biggr ], \end{align}

where $\zeta$ is defined in Eq. (5.41). In this way, we have reduced the problem of estimating $1/(\varepsilon _1 \sigma _0)$ to a simple numerical quadrature of the function $S(\zeta )$ in Eq. (5.42). To obtain a more explicit, but less accurate, approximation, we neglect the $a \pi ^2\xi ^2/9$ term in Eq. (5.42), which corresponds to using the result (B.12) for $\kappa _j$ , and then evaluate the resulting integral to get $S= {\left (\zeta - \ln (1+\zeta )\right )/\zeta ^2}$ . In this way, Eq. (5.38) can be approximated more explicitly as

(5.44) \begin{equation} \frac {1}{\varepsilon _1 \sigma _0} \approx \frac {2}{\pi } \left [\frac {\zeta ^2}{\zeta - \ln (1+\zeta )} + C_1 - \ln \left (\frac {2\varepsilon _1}{\varepsilon } \right )\right ]. \end{equation}

For $63$ Dirichlet patches (i.e. $N=64$ ) and with $\varepsilon _1=0.1$ , Figure 10 compares the asymptotic results obtained by using the discrete sum (5.38) with its large- $N$ approximation (5.43) and with the simpler, more explicit, result (5.44). We observe that Eq. (5.43) provides an excellent approximation, while Eq. (5.44) has a small systematic underestimate.

6.1. Interior targets

Throughout this paper, we focused on reactive patches on the boundary of a bounded domain. In many applications, however, absorbing sinks, traps and/or reactive targets can be hidden inside a bounded domain $\Omega _0$ , surrounded by a reflecting boundary $\partial \Omega _0$ . Let us consider the problem with $N$ interior targets, where each target is a compact set $\Omega _{\varepsilon _j} \subset \Omega _0$ of size $\varepsilon _j$ , centred at a point ${{\boldsymbol{x}}}_j$ (Figure 11). Since the targets are impenetrable for a diffusing particle, we still consider surface reactions on their boundaries, denoted as $\Gamma _{\varepsilon _j} = \partial \Omega _{\varepsilon _j}$ . As before, the targets are small ( $\varepsilon _j \sim o(1)$ ), comparable in size, and well-separated from each other: $|{{\boldsymbol{x}}}_j - {{\boldsymbol{x}}}_k| \sim {\mathcal O}(1)$ for $j\ne k$ , and from the domain boundary $\partial \Omega _0$ : $|{{\boldsymbol{x}}}_j - {{\boldsymbol{x}}}| \sim {\mathcal O}(1)$ for any ${{\boldsymbol{x}}}\in \partial \Omega _0$ . This setting represents diffusion in a perforated domain $\Omega = \Omega _0 \backslash (\Omega _{\varepsilon _1} \cup \cdots \cup \Omega _{\varepsilon _N})$ with the boundary $\partial \Omega = \partial \Omega _0 \cup \Gamma _{\varepsilon _1} \cup \cdots \cup \Gamma _{\varepsilon _N}$ . These notations allow us to make an equivalence with the earlier setting introduced at the beginning of Section 2. In particular, we can retain the same formulations of the four considered problems. As expected, their solutions can be constructed analogously, but with some modifications in the ‘building blocks’. In this section, we briefly describe these modifications and illustrate a few results.

Figure 11. Illustration of a bounded domain $\Omega \subset {\mathbb{R}}^2$ with a smooth reflecting boundary $\partial \Omega _0$ (grey dashed line) and three interior targets $\Omega _{\varepsilon _j}$ (filled in grey), centred at ${{\boldsymbol{x}}}_j$ , with reactive boundaries $\Gamma _{\varepsilon _j}$ (in red and blue). For a particle started from a point ${{\boldsymbol{x}}}\in \Omega$ , the splitting probability $S_1({{\boldsymbol{x}}})$ is the probability of hitting the boundary $\Gamma _{\varepsilon _1}$ first.

6.1.1 Perfectly reactive targets: Dirichlet boundary condition

If the $j$ -th target is perfectly reactive (with Dirichlet condition), the inner solution around this target is proportional to the exterior Dirichlet Green’s function

(6.1) \begin{equation} \Delta g_\infty = 0 \quad \textrm {in}\, {\mathbb{R}}^2 \backslash \Omega _{j}; \qquad g_\infty = 0 \quad \textrm {on}\, \partial \Omega _{j} ; \qquad g_\infty \sim \ln |{\boldsymbol{y}}| + {\mathcal O}(1) \quad \textrm {as}\, |{\boldsymbol{y}}| \to \infty , \end{equation}

where $\Omega _j = \varepsilon _j^{-1} \Omega _{\varepsilon _j}$ is the rescaled target $\Omega _{\varepsilon _j}$ . In contrast to the earlier studied case of a Dirichlet patch, the Green’s function $g_\infty (\,{\boldsymbol{y}})$ depends on the shape of the target $\Omega _j$ and thus is not universal. In particular, its asymptotic behaviour is $g_\infty (\,{\boldsymbol{y}}) \sim \ln |{\boldsymbol{y}}| - \ln (d_j) + o(1)$ as $|{\boldsymbol{y}}|\to \infty$ , where $d_j$ is the logarithmic capacity of $\Omega _j$ . For instance, if $\Omega _j$ is the unit disk, its logarithmic capacity is $1$ . Numerical values for $d_j$ for various shapes of $\Omega _j$ are given in Table 1 of [Reference Kurella, Tzou, Coombs and Ward54].

In addition, the outer solution involves the Neumann Green’s function $G_{\textrm {b}}({{\boldsymbol{x}}},{\boldsymbol \xi })$ (also known as pseudo-Green’s function) that satisfies

(6.2a) \begin{align} \Delta G_{\textrm {b}} & = \frac {1}{|\Omega _0|}-\delta ({{\boldsymbol{x}}}-{\boldsymbol \xi }) \quad \textrm {in}\, \Omega _0; \qquad \partial _n G_{\textrm {b}} = 0 \quad \textrm {on}\, \partial \Omega _0 ; \qquad \int \limits _{\Omega _0} G_{\textrm {b}}({{\boldsymbol{x}}},{\boldsymbol \xi })\, d{{\boldsymbol{x}}} = 0 , \\[-16pt] \nonumber \end{align}

(6.2b) \begin{align} G_{\textrm {b}}({{\boldsymbol{x}}},{\boldsymbol \xi }) & \sim - \frac {1}{2\pi } \ln |{{\boldsymbol{x}}}-{\boldsymbol \xi }| + R_{\textrm {b}}({\boldsymbol \xi }) + o(1) \quad \textrm {as}\, {{\boldsymbol{x}}} \to {\boldsymbol \xi }\in \Omega _0, \\[8pt] \nonumber \end{align}

where we added the subscript $b$ to distinguish it from the surface Neumann Green’s function $G({{\boldsymbol{x}}},{\boldsymbol \xi })$ . The main difference between this Green’s function and the surface Neumann Green’s function defined by Eqs. (2.9) is that the singularity at $\boldsymbol \xi$ is located in the bulk and not on the boundary (accordingly, there is the factor $1/(2\pi )$ in Eq. (6.2b) instead of $1/\pi$ ). For instance, when $\Omega _0$ is the unit disk, the Neumann Green’s function is well-known [Reference Kolokolnikov, Titcombe and Ward51]:

(6.3a) \begin{align} G_{\textrm {b}}({{\boldsymbol{x}}},{\boldsymbol \xi }) & = -\frac {1}{2\pi }\ln |{{\boldsymbol{x}}}-{\boldsymbol \xi }| - \frac {1}{4\pi } \ln \left (|{{\boldsymbol{x}}}|^2|{\boldsymbol \xi }|^2 + 1 -2{{\boldsymbol{x}}} \boldsymbol{\cdot }{\boldsymbol \xi } \right ) + \frac {|{{\boldsymbol{x}}}|^2 + |{\boldsymbol \xi }|^2}{4\pi } - \frac {3}{8\pi }, \\[-10pt] \nonumber \end{align}

(6.3b) \begin{align} R_{\textrm {b}}({\boldsymbol \xi }) & = -\frac {1}{2\pi }\ln \left (1-|{\boldsymbol \xi }|^2\right ) + \frac {|{\boldsymbol \xi }|^2}{2\pi } - \frac {3}{8\pi } . \\[8pt] \nonumber \end{align}

We remark that rapidly converging infinite series representations for $G_{\textrm {b}}$ and $R_{\textrm {b}}$ are also known explicitly for ellipses (see Eqs. (4.6, 4.7) from [Reference Iyaniwura, Wong, Macdonald and Ward50]) and for rectangles (see Eq. (4.13) of Section 4.2 of [Reference Kolokolnikov, Ward and Wei52], as well as [Reference McCann, Hazlett and Babu61]).

With these minor changes in the ‘building blocks’, we can repeat the steps in Sections 2.2 and 2.3 to obtain that the splitting probability $S_k({{\boldsymbol{x}}})$ is

(6.4) \begin{equation} S_k({{\boldsymbol{x}}}) = \chi _k - \sum \limits _{i=1}^N 2\pi A_i G_{\textrm {b}}({{\boldsymbol{x}}},{{\boldsymbol{x}}}_i) , \end{equation}

where the coefficients $\chi _k$ and $A_i$ are still determined via Eqs. (2.17, 2.21), but now with $\nu _j = -1/\ln (\varepsilon _jd_j)$ . Note also that the matrix $\textbf{G}$ from Eq. (2.15) is now based on the Neumann Green’s function $G_{\textrm {b}}({{\boldsymbol{x}}},{\boldsymbol \xi })$ from Eq. (6.2) and its regular part $R_{\textrm {b}}({\boldsymbol \xi })$ . In addition, the factor $\pi$ in the definition (2.14) of the elements of the matrix $\textbf{G}$ should be replaced by $2\pi$ . The example of two targets can be worked out explicitly. Moreover, equally spaced targets located on a circular ring that is concentric within a unit disk can also be treated explicitly.

6.1.2 Partially reactive targets: Robin boundary condition

If the $j$ -th target is partially reactive, we should replace the Dirichlet boundary condition by a Robin condition with the reactivity parameter $q_j$ . In the same vein, the former Robin Green’s function $g_\mu (\,{\boldsymbol{y}})$ now satisfies

(6.5) \begin{equation} \Delta g_{\mu } = 0 \quad \textrm {in}\, {\mathbb{R}}^2 \backslash \Omega _{j}; \qquad \partial _n g_{\mu } + \mu g_{\mu } = 0 \quad \textrm {on}\, \partial \Omega _{j} ; \qquad g_{\mu } \sim \ln |{\boldsymbol{y}}| + {\mathcal{C}}_j(\varepsilon _j q_j) + o(1) \quad \textrm {as}\, |{\boldsymbol{y}}| \to \infty , \end{equation}

with $\mu = \varepsilon _j q_j$ . In particular, the constant term ${\mathcal{C}}_j(\mu )$ of its asymptotic behaviour at infinity is not universal and depends on the shape of $\Omega _j$ . We can still employ the eigenmodes of the exterior Steklov problem in ${\mathbb{R}}^2 \backslash \Omega _j$ to construct $g_\mu (\,{\boldsymbol{y}})$ and to determine the spectral expansion for ${\mathcal{C}}_j(\mu )$ :

(6.6) \begin{equation} {\mathcal{C}}_j(\mu ) = -\ln (d_j) + \frac {2\pi }{\mu |\partial \Omega _j|} + 2\pi \sum \limits _{k=1}^\infty \frac {[\Psi _k^{\,j}(\infty )]^2}{\mu _k^{\,j} + \mu } , \end{equation}

where $\mu _k^{\,j}$ and $\Psi _k^{\,j}$ are the eigenvalues and eigenfunctions of the auxiliary exterior Steklov problem,

(6.7) \begin{equation} \Delta \Psi _k^j = 0 \quad \textrm {in}\, {\mathbb{R}}^2 \backslash \Omega _j; \qquad \partial _n \Psi _k^j = \mu _k^j \Psi _k^j \quad \textrm {on}\, \partial \Omega _j ; \qquad \Psi _k^j \sim {\mathcal O}(1) \quad \textrm {as}\, |{\boldsymbol{y}}| \to \infty . \end{equation}

We also used the normalisation $\Psi _0^j = 1/\sqrt {|\partial \Omega _j|}$ of the principal eigenfunction associated to $\mu _0^j = 0$ . As earlier, a partially reactive target can be treated a perfect one, but with the reduced size:

(6.8) \begin{equation} \varepsilon _j^{\textrm {eff}} = \varepsilon _j \,\exp \bigl ({-}\ln (d_j) - {\mathcal{C}}_j(\varepsilon _j q_j)\bigr ), \end{equation}

(see further discussion and examples in Section 3).

For the target $\Omega _j$ of an arbitrary shape, the computation of the Steklov eigenmodes $\mu _k^j$ and $\Psi _k^j$ requires numerical techniques (e.g., a finite-element method, see [Reference Chaigneau and Grebenkov11, Reference Grebenkov and Chaigneau39] and references therein). However, if the rescaled target $\Omega _j$ is the unit disk, the eigenmodes are known explicitly and they all vanish at infinity, except $\Psi _0^j$ . As the logarithmic capacity of the unit disk is equal to $1$ , we get a particularly simple exact expression:

(6.9) \begin{equation} {\mathcal{C}}_{\textrm {disk}}(\mu ) = \frac {1}{\mu } . \end{equation}

This is not surprising given that the exterior Green’s function for the unit disk is simply $g_\mu (\,{\boldsymbol{y}}) = 1/\mu + \ln |{\boldsymbol{y}}|$ .

6.1.3 Imperfect targets

In a similar way, we can easily reproduce the derivations of Sections 4 and 5 for the mixed Steklov-Neumann and Steklov-Neumann-Dirichlet problems. For an imperfect target $\Omega _{\varepsilon _j}$ of arbitrary shape, the main difficulty is the lack of knowledge of the function ${\mathcal{C}}_j(\mu )$ , which is formally accessible via the spectral expansion (6.6) but its ‘ingredients’ require numerical computations. Moreover, as the coefficients of the Taylor expansion of ${\mathcal{C}}_j(\mu )$ are unknown, we cannot rely on the approximation (3.10). As a consequence, many numerical steps would be involved, and the analytical, almost explicit form of the asymptotic results would in general be lost. An interesting extension of this work consists of a systematic study of the function ${\mathcal{C}}_j(\mu )$ for targets of various shapes.

A drastic simplification appears when the imperfect targets are disks due the explicit form (6.9) of the function ${\mathcal{C}}(\mu )$ . In this case, the analysis of Sections 4 and 5 can be reproduced and will actually be even simpler. For instance, for the mixed Steklov-Neumann-Dirichlet problem with a single disk-shaped target $\Omega _{\varepsilon _1}$ and one perfectly reactive target $\Omega _{\varepsilon _2}$ (of arbitrary shape), one can rewrite Eq. (5.28) as

(6.10) \begin{equation} C = \ln (d_2 \varepsilon _1 \varepsilon _2) - 2\pi \bigl [R_{\textrm {b}}({{\boldsymbol{x}}}_1) + R_{\textrm {b}}({{\boldsymbol{x}}}_2) - 2G_{\textrm {b}}({{\boldsymbol{x}}}_1,{{\boldsymbol{x}}}_2)\bigr ], \end{equation}

from which the principal eigenvalue reads

(6.11) \begin{equation} \frac {1}{\varepsilon _1 \sigma _0} \approx - \ln (d_2 \varepsilon _1 \varepsilon _2) + 2\pi \bigl [R_{\textrm {b}}({{\boldsymbol{x}}}_1) + R_{\textrm { b}}({{\boldsymbol{x}}}_2) - 2G_{\textrm {b}}({{\boldsymbol{x}}}_1,{{\boldsymbol{x}}}_2)\bigr ]. \end{equation}

Figure 12 illustrates the accuracy of the asymptotic relation (6.11) for the unit disk with two interior circular targets. We observe a close agreement between a numerical solution and the asymptotic formula. One can notice a small deviation between two lines that slightly increases as the target radius $\varepsilon _2$ decreases. This minor discrepancy seems to be a numerical artefact due to the available mesh size $0.005$ , which becomes comparable to the target at small $\varepsilon _2$ . To check this point, we computed the eigenvalue of the mixed Steklov-Dirichlet problem for a circular annulus with radii $\varepsilon _2$ and $R = 1$ , with Steklov condition on the outer circle and Dirichlet condition on the inner circle. As the exact solution of this problem is known, $\sigma _0 = 1/(R\ln (R/\varepsilon _2))$ , we could compare it with the numerical results, and found the same minor discrepancy.

Figure 12. The asymptotic behaviour of the principal eigenvalue of the mixed Steklov- Neumann-Dirichlet problem, plotted as $1/(\varepsilon _1\sigma _0)$ versus $\varepsilon _2$ , for the unit disk with two interior circular targets of radii $\varepsilon _1 = 0.05$ and $\varepsilon _2$ (variable from $0.01$ to $0.1$ ), located at ${{\boldsymbol{x}}}_1 = ({-}0.5,0)$ and ${{\boldsymbol{x}}}_2 = (0.5,0)$ . Filled circles present the numerical solution by a FEM with the maximal mesh size of $0.005$ , solid line indicates Eq. (6.11).

In summary, we conclude that interior targets can be handled in essentially the same way as boundary patches, even though the asymptotic analysis becomes sensitive to the shapes of the targets. Moreover, one can combine interior targets with boundary patches that opens a way to access a broad variety of various geometric settings.

6.2. Exterior problems

Another extension of the present approach is related to exterior problems in $\Omega = {\mathbb{R}}^2 \backslash \Omega _0$ , where $\Omega _0$ is a simply-connected compact domain. While the inner solutions remain unchanged, the outer solution is now constructed using the exterior surface Neumann Green’s function, labelled by $G_{\textrm {e}}$ , which satisfies

(6.12a) \begin{align} \Delta G_{\textrm {e}} & = 0 \quad \textrm {in}\, {\mathbb{R}}^2\backslash \Omega _0, \\[-12pt] \nonumber \end{align}

(6.12b) \begin{align} G_{\textrm {e}}({{\boldsymbol{x}}},{\boldsymbol \xi }) & \sim - \frac {1}{\pi } \ln |{{\boldsymbol{x}}}-{\boldsymbol \xi }| + R_{\textrm {e}}({\boldsymbol \xi }) + o(1) \quad \textrm {as}\, {{\boldsymbol{x}}} \to {\boldsymbol \xi }\in \partial \Omega _0; \qquad \partial _n G_{\textrm {e}} = 0 \quad \textrm {on}\, \partial \Omega _0 \backslash \{{\boldsymbol \xi }\} , \\[-12pt] \nonumber \end{align}

(6.12c) \begin{align} G_{\textrm {e}} ({{\boldsymbol{x}}},{\boldsymbol \xi }) & \sim -\frac {1}{2\pi }\ln |{{\boldsymbol{x}}}| + o(1) \quad \textrm {as}\, |{{\boldsymbol{x}}}|\to \infty , \\[8pt] \nonumber \end{align}

where $R_{\textrm {e}}({\boldsymbol \xi })$ is the regular part of $G_{\textrm {e}}({{\boldsymbol{x}}},{\boldsymbol \xi })$ at $\boldsymbol \xi$ . The condition that $G_{\textrm {e}}({{\boldsymbol{x}}},{\boldsymbol \xi }) +(2\pi )^{-1}\ln |{{\boldsymbol{x}}}|\to 0$ as $|{{\boldsymbol{x}}}|\to \infty$ determines $G_{\textrm {e}}$ uniquely.

In the case when $\Omega _0$ is the unit disk and $|{\boldsymbol \xi }|=1$ , we claim that

(6.13) \begin{equation} G_{\textrm {e}}({{\boldsymbol{x}}},{\boldsymbol \xi }) = -\frac {1}{\pi }\ln |{{\boldsymbol{x}}}-{\boldsymbol \xi }| + \frac {1}{2\pi } \ln |{{\boldsymbol{x}}}| , \qquad R_{\textrm {e}}({\boldsymbol \xi }) = 0 . \end{equation}

Clearly Eq. (6.13) satisfies the behaviour (6.12c) as $|{{\boldsymbol{x}}}|\to \infty$ as well as Eq. (6.12b) as ${{\boldsymbol{x}}}\to {\boldsymbol \xi }$ , where we identify that $R_{\textrm {e}}({\boldsymbol \xi }) = 0$ (see more details in Appendix F.1).

For instance, if there are two patches on the unit circle, the solution of any of four earlier considered problems involves

(6.14) \begin{equation} R_{\textrm {e}}({{\boldsymbol{x}}}_1) + R_{\textrm {e}}({{\boldsymbol{x}}}_2) - 2G_{\textrm {e}}({{\boldsymbol{x}}}_1,{{\boldsymbol{x}}}_2) = \frac {2}{\pi } \ln |{{\boldsymbol{x}}}_1-{{\boldsymbol{x}}}_2|, \end{equation}

which is identical for both interior and exterior domains. This property is consistent with the fact that the eigenvalues of interior and exterior mixed Steklov problems for the unit disk are identical. In contrast, the associated eigenfunctions behave differently.