# Spectral properties of the Dirichlet-to-Neumann operator for spheroids

Denis S Grebenkov

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

## To cite this version:

Denis S Grebenkov. Spectral properties of the Dirichlet-to-Neumann operator for spheroids. Physical Review E, 2024, 109 (5), pp.055306. 10.1103/PhysRevE.109.055306. hal-04749302

HAL Id: hal-04749302

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

Submitted on 23 Oct 2024

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

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

<figure><img src="image_2.png" /></figure>
Spectral properties of the Dirichlet-to-Neumann operator for spheroids

Denis S. Grebenkov¹, *

¹Laboratoire de Physique de la Matière Condensée (UMR 7643),
CNRS - Ecole Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

(Dated: April 16, 2024)

We study the spectral properties of the Dirichlet-to-Neumann operator and the related Steklov problem in spheroidal domains ranging from a needle to a disk. An explicit matrix representation of this operator for both interior and exterior problems is derived. We show how the anisotropy of spheroids affects the eigenvalues and eigenfunctions of the operator. As examples of physical applications, we discuss diffusion-controlled reactions on spheroidal partially reactive targets and the statistics of encounters between the diffusing particle and the spheroidal boundary.

PACS numbers: 02.50.-r, 05.40.-a, 02.70.Rr, 05.10.Gg

Keywords: diffusion, spheroids, Steklov problem, Dirichlet-to-Neumann operator, Laplace equation

# I. INTRODUCTION

The Dirichlet-to-Neumann operator $M$ plays an important role in applied mathematics, physics, and engineering. One of its most known applications is related to electrical impedance tomography, also known as Calderon's problem [1-5], in which the electric conductivity in the bulk has to be determined from electric measurements on the boundary. The Dirichlet-to-Neumann operator is also employed as a “building block” for analyzing and solving spectral and scattering problems in complex media via domain decomposition (see [6-10] and references therein). The eigenfunctions of $M$ often appear as a basis for representing and approximating harmonic functions and related quantities [11-15]. In chemical physics, diffusion-controlled reactions and other diffusion-mediated surface phenomena can be described by means of the operator $M$ [16-19]. In particular, the statistics of encounters between a diffusing particle and the confining boundary can be determined via a spectral expansion based on the Dirichlet-to-Neumann operator [20-22]. Its relations to other first-passage time statistics were also investigated [18, 23].

The spectral properties of the Dirichlet-to-Neumann operator and the related Steklov problem were thoroughly investigated for Euclidean domains and Riemannian manifolds (see the recent book [24] and reviews [25, 26]). We focus here on the most basic setting of an Euclidean domain $\Omega \subset \mathbb{R}^d$ with a smooth bounded boundary $\partial\Omega$. The Dirichlet-to-Neumann operator $M$ is defined as a map of a function $f \in H^{\frac{1}{2}}(\partial\Omega)$ on the boundary $\partial\Omega$ onto another function $g \in H^{-\frac{1}{2}}(\partial\Omega)$ on that boundary such that $g = (\partial_n u)|_{\partial\Omega} = Mf$, where $\partial_n$ is the normal derivative oriented outwards the domain $\Omega$, and $u$ is the unique solution of the Dirichlet boundary value problem:

$$ \Delta u = 0 \quad \text{in } \Omega, \quad u = f \quad \text{on } \partial\Omega, \quad (1) $$

with $\Delta$ being the Laplace operator (here $H^{\pm \frac{1}{2}}(\partial\Omega)$ are appropriate functional spaces, see Ref. [24] for mathematical details and references). In other words, the operator $M$ transforms the Dirichlet boundary condition $u|_{\partial\Omega} = f$ on $\partial\Omega$ into an equivalent Neumann boundary condition $(\partial_n u)|_{\partial\Omega} = Mf$. For instance, if $f$ is a given concentration or temperature profile maintained on $\partial\Omega$, then the Laplace equation describes the steady-state regime of molecular or heat diffusion, and $(\partial_n u)|_{\partial\Omega}$ is proportional to the flux density on $\partial\Omega$.

When $\Omega$ is bounded, $M$ is known to be pseudo-differential self-adjoint operator with a discrete spectrum, i.e., there is an infinite countable sequence of eigen-pairs $\{\mu_k, v_k\}$ satisfying $Mv_k = \mu_k v_k$; the nonnegative eigenvalues $\mu_k$ are enumerated by $k = 0, 1, 2, \dots$ in an increasing order,

$$ 0 \le \mu_0 \le \mu_1 \le \dots \nearrow +\infty, \quad (2) $$

whereas the associated eigenfunctions $\{v_k\}$ form a complete basis in the space $L^2(\partial\Omega)$ of square-integrable functions on $\partial\Omega$ [24]. Alternatively, one can search for solutions of the Steklov problem,

$$ \Delta V_k = 0 \quad \text{in } \Omega, \quad \partial_n V_k = \mu_k V_k \quad \text{on } \partial\Omega, \quad (3) $$

where the Steklov eigenvalues standing in the boundary condition are identical to $\mu_k$. This tight relation implies that each Steklov eigenfunction $V_k$ can be obtained as a harmonic extension of the eigenfunction $v_k$ of $M$.

Despite numerous mathematical studies of spectral properties of the Dirichlet-to-Neumann operator [24], intricate relations between its spectrum and the geometric features of the boundary $\partial\Omega$ are not yet fully understood. The eigenvalues and eigenfunctions of $M$ are known explicitly only in few simple domains such as a ball, a space between concentric spheres, the exterior of a ball, and rectangular cuboids [24, 27]. In particular, the role of the boundary anisotropy remains unclear. The situation is even worse for the exterior problem when $\Omega = \mathbb{R}^d \setminus \Omega_0$ is the exterior of a bounded domain $\Omega_0$. Even though the domain $\Omega$ is unbounded, its boundary $\partial\Omega$ is bounded that implies again the discrete spectrum of the Dirichlet-to-Neumann operator $M$ [12, 13, 15, 28-30]. However,

<footnote>* denis.grebenkov@polytechnique.edu</footnote>
the analysis of the exterior problem is more difficult; in particular, the mathematical proofs substantially differ for space dimensions $d = 2$ and $d \ge 3$. Relations between spectral properties and the geometric shape of $\partial\Omega$ were much less studied. For instance, to our knowledge, the exterior of a ball is the unique example, for which the eigenvalues and eigenfunctions of $M$ are known explicitly for the exterior problem.

In this paper, we study the spectral properties of the Dirichlet-to-Neumann operator $M$ on prolate and oblate spheroidal surfaces that allow one to model various anisotropic shapes in three dimensions, ranging from a needle to a disk. We focus on the less studied exterior spectral problem (an extension to the interior problem is summarized in Appendix A). By employing the prolate/oblate spheroidal coordinates to represent a general solution of the Laplace equation (1), we obtain a convenient matrix representation of the operator $M$. This matrix can then be truncated and diagonalized numerically to approximate the eigenvalues $\mu_k$ and eigenfunctions $v_k$ of $M$. This efficient technique allows us to investigate how the spectral properties of the Dirichlet-to-Neumann operator depend on the anisotropy of the boundary, especially in the limits of elongated (needle-like) and flattened (disk-like) spheroids. While similar techniques were applied in the past for solving various boundary value problems in spheroidal domains (see, e.g., [31–42] and references therein), we are not aware of earlier studies of the spectral properties of the Dirichlet-to-Neumann operator in these domains. Sections II and III are devoted respectively to prolate and oblate spheroidal domains. In Sec. IV, we discuss two applications of these results for understanding diffusion-controlled reactions and the statistics of boundary encounters. Section V summarizes our findings and presents future perspectives.

## II. PROLATE SPHEROIDS

In this section, we study the Dirichlet-to-Neumann operator $M$ in the exterior of a prolate spheroid with semi-axes $a \le b$:

$$ \Omega = \left\{ (x, y, z) \in \mathbb{R}^3 : \frac{x^2}{a^2} + \frac{y^2}{a^2} + \frac{z^2}{b^2} > 1 \right\}. \quad (4) $$

In the prolate spheroidal coordinates $(\alpha, \theta, \phi)$,

$$ \begin{pmatrix} x \\ y \\ z \end{pmatrix} = a_E \begin{pmatrix} \sinh \alpha \sin \theta \cos \phi \\ \sinh \alpha \sin \theta \sin \phi \\ \cosh \alpha \cos \theta \end{pmatrix} \quad \begin{cases} 0 < \alpha < \infty \\ 0 \le \theta \le \pi \\ 0 \le \phi < 2\pi \end{cases}, $$

with $a_E = \sqrt{b^2 - a^2}$ (Fig. 1a), the scale factors determining the surface and volume elements, are [43]

$$ h_{\alpha} = h_{\theta} = a_E \sqrt{\sinh^2 \alpha + \sin^2 \theta}, \quad (5a) $$

$$ h_{\phi} = a_E \sinh \alpha \sin \theta. \quad (5b) $$

<figure><img src="image_3.png" /><figcaption>FIG. 1. Illustration for prolate (a) and oblate (b) spheroidal coordinates $(\alpha, \theta, \phi)$. Note that the angle $\theta$ is defined differently in two cases.</figcaption></figure>

In these coordinates, the domain $\Omega$ is characterized by $\alpha > \alpha_0 = \tanh^{-1}(a/b)$, while its boundary $\partial\Omega$ is determined by the condition $\alpha = \alpha_0$. The action of the Laplace operator onto a function $u$ reads

$$ \Delta u = \frac{1}{a_E^2(\sinh^2 \alpha + \sin^2 \theta)} \left[ \frac{1}{\sinh \alpha} \frac{\partial}{\partial \alpha} \left( \sinh \alpha \frac{\partial u}{\partial \alpha} \right) + \frac{1}{\sin \theta} \frac{\partial}{\partial \theta} \left( \sin \theta \frac{\partial u}{\partial \theta} \right) \right] + \frac{1}{a_E^2 \sinh^2 \alpha \sin^2 \theta} \frac{\partial^2 u}{\partial \phi^2}. \quad (6) $$

### Matrix representation

A general solution of the Laplace equation reads

$$ u(\alpha, \theta, \phi) = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} A_{mn} Q_n^m(\cosh \alpha) Y_{mn}(\theta, \phi), \quad (7) $$

where $A_{mn}$ are arbitrary coefficients, $Q_n^m(z)$ are the associated Legendre functions of the second kind (see Appendix B), and

$$ Y_{mn}(\theta, \phi) = a_{mn} P_n^m(\cos \theta) e^{im\phi}, \quad (8) \\ a_{mn} = \sqrt{\frac{2n+1}{4\pi} \frac{(n-m)!}{(n+m)!}} $$

are the normalized spherical harmonics, with $P_n^m(x)$ being the associated Legendre polynomials. The normal derivative reads

$$ \frac{\partial u}{\partial n} \Big|_{\partial \Omega} = - \frac{1}{h_{\alpha}} \partial_{\alpha} u \Big|_{\alpha=\alpha_0} \quad (9) \\ = - \sum_{n=0}^{\infty} \sum_{m=-n}^{n} A_{mn} Y_{mn}(\theta, \phi) \frac{\sinh \alpha_0 Q_n^{m'}(\cosh \alpha_0)}{a_E \sqrt{\cosh^2 \alpha_0 - \cos^2 \theta}}, $$

where prime denotes the derivative with respect to the argument. In other words, for a square-integrable function $f$ on $\partial\Omega$, decomposed on the complete basis of spherical harmonics,

$$ f(\theta, \phi) = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} f_{mn} Y_{mn}(\theta, \phi) \quad (10) $$
(with coefficients $f_{mn} = A_{mn}Q_n^m(\cosh \alpha_0)$), the action of the Dirichlet-to-Neumann operator $\mathcal{M}$ is

$$[\mathcal{M}f](\theta, \phi) = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} f_{mn} \frac{c_{mn} Y_{mn}(\theta, \phi)}{\sqrt{\cosh^2 \alpha_0 - \cos^2 \theta}}, \quad (11)$$

where

$$c_{mn} = - \frac{\sinh \alpha_0 Q_n^m (\cosh \alpha_0)}{a_E Q_n^m (\cosh \alpha_0)}. \quad (12)$$

At the same time, the function $\mathcal{M}f$ can also be decomposed on the complete basis of $Y_{mn}$:

$$[\mathcal{M}f](\theta, \phi) = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} g_{mn} Y_{mn}(\theta, \phi). \quad (13)$$

Its coefficients $g_{mn}$ can be obtained by multiplying Eq. (11) by $Y_{m'n'}^*(\theta, \phi) \sin \theta$ and integrating over $\theta$ and $\phi$:

$$\begin{aligned} g_{m'n'} &= \int_0^{2\pi} d\phi \int_0^{\pi} d\theta \sin\theta Y_{m'n'}^*(\theta, \phi) \\ &\quad \times \sum_{n=0}^{\infty} \sum_{m=-n}^{n} f_{mn} \frac{c_{mn} Y_{mn}(\theta, \phi)}{\sqrt{\cosh^2 \alpha_0 - \cos^2 \theta}} \\ &= \sum_{n=0}^{\infty} \sum_{m=-n}^{n} f_{mn} \mathbf{M}_{mn,m'n'} \end{aligned}$$

where

$$\mathbf{M}_{mn,m'n'} = 2\pi \delta_{m,m'} a_{mn} a_{mn'} c_{m'n'} F_{n,n'}^m (\cosh \alpha_0), \quad (14)$$

and

$$F_{n,n'}^m(z) = \int_{-1}^{1} dx \frac{P_n^m(x) P_{n'}^m(x)}{\sqrt{z^2 - x^2}}. \quad (15)$$

We stress that the elements of the (infinite-dimensional) matrix $\mathbf{M}$ are enumerated by double indices $mn$ and $m'n'$. In practice, we employ the following order:

$$mn: 00 (-1)1 01 11 (-2)2 (-1)2 02 12 22 \dots$$

which is borrowed from the enumeration of spherical harmonics. As a consequence, Eq. (13) can be written as

$$[\mathcal{M}f](\theta, \phi) = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} [\mathbf{M}f]_{mn} Y_{mn}(\theta, \phi), \quad (16)$$

i.e., the matrix $\mathbf{M}$ represents the operator $\mathcal{M}$ on the orthonormal basis of spherical harmonics.

The eigenvalues of the matrix $\mathbf{M}$ coincide with the eigenvalues $\mu_k$ of the Dirichlet-to-Neumann operator $\mathcal{M}$; in turn, each eigenvector $\mathbf{V}_k$ of $\mathbf{M}$, satisfying $\mathbf{M}\mathbf{V}_k = \mu_k \mathbf{V}_k$, determines the coefficients in the representation of

the associated eigenfunction $v_k$ on the basis of spherical harmonics:

$$v_k(\theta, \phi) = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} [\mathbf{V}_k]_{mn} Y_{mn}(\theta, \phi). \quad (17)$$

In Appendix B 4, we check the orthogonality of these eigenfunctions to each other; moreover, they can also be normalized as

$$\begin{aligned} (v_k, v_{k'})_{L^2(\partial\Omega)} &= \int_{\partial\Omega} ds v_k^* v_{k'} \\ &= \int_0^\pi d\theta \int_0^{2\pi} d\phi h_\theta h_\phi v_k^*(\theta, \phi) v_{k'}(\theta, \phi) = \delta_{k,k'}. \end{aligned} \quad (18)$$

Using Eq. (7) for a general solution of the Laplace operator, one easily finds a harmonic extension of the eigenfunction $v_k$ into $\Omega$, i.e., the Steklov eigenfunction associated to the eigenvalue $\mu_k$:

$$V_k(\alpha, \theta, \phi) = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} \frac{Q_n^m(\cosh \alpha)}{Q_n^m(\cosh \alpha_0)} [\mathbf{V}_k]_{mn} Y_{mn}(\theta, \phi). \quad (19)$$

The asymptotic behavior of $Q_n^m(z)$ for large $|z|$, $Q_n^m(z) \propto z^{-n-1}$ implies the expected power-law decay of Steklov eigenfunctions in the leading order as $|x| \to \infty$:

$$\begin{aligned} V_k(\alpha, \theta, \phi) &\approx \frac{[\mathbf{V}_k]_{00}}{\sqrt{4\pi} Q_0(\cosh \alpha_0)} \frac{1}{\cosh \alpha} \\ &\approx \frac{[\mathbf{V}_k]_{00}}{\sqrt{4\pi} Q_0(\cosh \alpha_0)} \frac{a_E}{|x|} \propto |x|^{-1}, \end{aligned} \quad (20)$$

where $|x| = a_E \sqrt{\cosh^2 \alpha - \sin^2 \theta}$ is the distance from the origin to a point $x = (\alpha, \theta, \phi)$.

## Classification of eigenfunctions

Since the Dirichlet-to-Neumann operator $\mathcal{M}$ does not affect the angle $\phi$, the matrix elements $\mathbf{M}_{mn,m'n'}$ are nonzero only when $m = m'$. In other words, the action of $\mathcal{M}$ onto a function $f(\theta)e^{im\phi}$ does not alter its dependence on $\phi$: $\mathcal{M}(f(\theta)e^{im\phi}) = g(\theta)e^{im\phi}$. As a consequence, any eigenfunction $v_k$ depends on $\phi$ via a factor $e^{im\phi}$ for some integer $m$ (or via a linear combination of $e^{im\phi}$ and $e^{-im\phi}$, see below). This property allows one to classify all eigenfunctions according to their dependence on $\phi$ and thus to enumerate them as $v_{mn}(\theta, \phi)$, in analogy to spherical harmonics $Y_{mn}(\theta, \phi)$. Here the index $m$ determines the dependence of the eigenfunction on $\phi$, $v_{mn}(\theta, \phi) \propto e^{im\phi}$, whereas the nonnegative index $n = |m|, |m| + 1, |m| + 2, \dots$ enumerates all such functions so that the associated eigenvalues $\mu_{mn}$ appear in an increasing order (for each fixed $m$):

$$0 \le \mu_{m|m|} \le \mu_{m(|m|+1)} \le \dots \quad (21)$$
Note that the index $n$ starts from $|m|$ in order to automatically satisfy the conventional restriction $|m| \le n$, known for spherical harmonics.

An alternative way to look at this classification consists in representing the matrix $\mathbf{M}$ as

$$ \mathbf{M} = \sum_{m=-\infty}^{\infty} \mathbf{M}_m, \qquad (22) $$

where the (sub)matrix $\mathbf{M}_m$ is composed of elements $\mathbf{M}_{mn,mn'}$ (and 0 otherwise). For instance, one has<fnref n="1" />

$$ \mathbf{M}_0 = \begin{pmatrix} \mathbf{M}_{00,00} & 0 & 0 & 0 & 0 & 0 & \mathbf{M}_{00,02} & 0 & 0 & \dots \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\ 0 & 0 & \mathbf{M}_{01,01} & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\ \mathbf{M}_{02,00} & 0 & 0 & 0 & 0 & 0 & \mathbf{M}_{02,02} & 0 & 0 & \dots \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \end{pmatrix}. \qquad (23) $$

If a vector $\mathbf{X}$ has the form $\mathbf{X}_{mn} = \delta_{m,m_0}x_n$ (i.e., its nonzero elements appear only for $m = m_0$), then $\mathbf{M}_m\mathbf{X} = 0$ for all $m \ne m_0$. In other words, the space $\ell^2$ of all vectors with square-summable elements can be decomposed into an (infinite) direct product of subspaces $\ell_m^2$ enumerated by $m$ ranging from $-\infty$ to $+\infty$. As a consequence, one can diagonalize separately each matrix $\mathbf{M}_m$ and then combine their eigenvalues $\mu_{mn}$ and eigenvectors $\mathbf{V}_{mn}$ (enumerated by the index $n = |m|, |m|+1, \dots)$ to construct the eigenvalues and eigenvectors of the matrix $\mathbf{M}$. We expect that the union of all eigenvalues $\mu_{mn}$ gives all eigenvalues of the matrix $\mathbf{M}$, i.e., such a decomposition determines the whole spectrum of $\mathbf{M}$. This statement is elementary in the finite-dimensional case, in particular, for a truncation of the matrix $\mathbf{M}$ that we will use for numerical computations. We can then rewrite the expansions (17, 19) as

$$ v_{mn}(\theta, \phi) = \sum_{n'=|m|}^{\infty} [\mathbf{V}_{mn}]_{mn'} Y_{mn'}(\theta, \phi), \qquad (24a) $$

$$ V_{mn}(\alpha, \theta, \phi) = \sum_{n'=|m|}^{\infty} \frac{Q_{n'}^m(\cosh \alpha)}{Q_{n'}^m(\cosh \alpha_0)} [\mathbf{V}_{mn}]_{mn'} Y_{mn'}(\theta, \phi). \qquad (24b) $$

We note that the symmetry

$$ P_n^m(-x) = (-1)^{m+n} P_n^m(x), \qquad (25) $$

implies $\mathbf{M}_{-m} = \mathbf{M}_m$ that allows one to restrict $m$ to be nonnegative. As a consequence, any eigenvalue $\mu_{mn}$ with

$m \ne 0$ should be (at least) twice degenerate. This degeneracy implies that any linear combination of $v_{(-m)n}$ and $v_{mn}$ is also an eigenfunction. As a consequence, if the eigenfunctions are constructed by using the decomposition (17) based on the diagonalization of the whole matrix $\mathbf{M}$ and then classified according to their dependence on $\phi$, each eigenfunction $v_k$ may in general exhibit the dependence on $\phi$ as a linear combination of $e^{im\phi}$ and $e^{-im\phi}$. An appropriate rotation by a $2 \times 2$ matrix can transform a pair of such eigenfunctions into those that are proportional to $e^{im\phi}$ and $e^{-im\phi}$. However, this step is not needed in practice, as we will diagonalize the matrices $\mathbf{M}_m$ to produce directly the desired dependence $e^{im\theta}$ on $\phi$.

In the following, we use interchangeably both notations $\mu_k, v_k$ and $\mu_{mn}, v_{mn}$ for eigenvalues and eigenfunctions. We recall that the single-index enumeration relies on the global ordering of all eigenvalues $\mu_k$ in Eq. (2). In turn, the double-index enumeration is based on the symmetries of eigenfunctions, namely, on their dependence on $\phi$ via $e^{im\phi}$, whereas the second index $n$ employs the ordering of $\mu_{mn}$ for each $m$ in Eq. (21).

## Limit of a sphere

In the limit $a \to b$, the prolate spheroid approaches the sphere of radius $b$. In this limit, one has $a_E \to 0$ and $\alpha_0 \to \infty$ such that $a_E \cosh \alpha_0 = b$ remains constant. As a consequence, the coefficients in Eq. (12) diverge as $c_{mn} \approx \frac{n+1}{b} \sinh \alpha_0$, whereas the matrix elements in Eq. (15) behave in the leading order as

$$ F_{n,n'}^{m}(\cosh \alpha_0) \approx \frac{\delta_{n,n'}}{2\pi a_{mn}^2 \cosh \alpha_0}, \qquad (26) $$

implying that $\mathbf{M}_{mn,m'n'} \to \frac{n+1}{b} \delta_{n,n'} \delta_{m,m'}$. The diagonal structure of this matrix yields

$$ \mu_{mn} = \frac{n+1}{b}, \qquad v_{mn} = \frac{1}{b} Y_{mn}(\theta, \phi). \qquad (27) $$

We retrieve therefore the well-known eigenvalues and eigenfunctions of the Dirichlet-to-Neumann operator for the exterior of a sphere. Note that the eigenvalues do not depend on $m$; moreover, since $m$ ranges from $-n$ to $n$, the degeneracy of the eigenvalue $\mu_{mn}$ is $2n+1$, as expected for a sphere due to its rotational symmetry. As illustrated below, the anisotropy of spheroids breaks this symmetry and reduces the degeneracy of eigenvalues.

## Numerical implementation

For a practical implementation, the infinite-dimensional matrix $\mathbf{M}$ has to be truncated to a finite size. In a basic setup, one can choose the truncation order $n_{max}$ to keep $n = 0, 1, 2, \dots, n_{max}$, and then construct the truncated matrix of size $(n_{max}+1)^2 \times (n_{max}+1)^2$,

<footnote n="1">¹ Note that $\mathbf{M}_{0n,0n'} = 0$ for odd $n + n'$ due to the symmetry (25) of the associated Legendre polynomials and the fact that the function $1/\sqrt{z^2 - x^2}$ in Eq. (15) is symmetric.</footnote>
as detailed in Appendix B. A numerical diagonalization of the truncated matrix provides an approximation for a number of eigenvalues and eigenfunctions of $\mathcal{M}$. As illustrated below, the accuracy of this approximation increases rapidly with the truncation order $n_{max}$.

A much faster procedure consists in dealing with the reduced matrices $\hat{\mathbf{M}}_m$, which are obtained from $\mathbf{M}_m$ by removing zero columns and rows. For instance, the matrix $\mathbf{M}_0$ from Eq. (23) has the following reduced form

$$ \hat{\mathbf{M}}_0 = \begin{pmatrix} \mathbf{M}_{00,00} & 0 & \mathbf{M}_{00,02} & \dots \\ 0 & \mathbf{M}_{01,01} & 0 & \dots \\ \mathbf{M}_{02,00} & 0 & \mathbf{M}_{02,02} & \dots \\ \dots & \dots & \dots & \dots \end{pmatrix}. \quad (28) $$

In practice, we start from the truncated matrix $\mathbf{M}$ of size $(n_{max}+1)^2 \times (n_{max}+1)^2$ and then dispatch its columns and rows according to $m$ ranging from 0 to $n_{max}$, into matrices $\hat{\mathbf{M}}_m$. As the reduced matrix $\hat{\mathbf{M}}_m$ is of much smaller size $(n_{max}+1-m) \times (n_{max}+1-m)$, its numerical diagonalization is significantly faster. Its eigenvalues approximate the eigenvalues $\mu_{mn}$ of the Dirichlet-to-Neumann operator $\mathcal{M}$; in turn, its eigenvectors determine the eigenfunctions $v_{mn}$ of $\mathcal{M}$ via Eq. (24a), and the Steklov eigenfunctions $V_{mn}$ via Eq. (24b). The former eigenfunctions are orthogonal to each other by construction. In turn, one needs to impose their normalization according to Eq. (18). We recall that this normalization is fixed up to an arbitrary phase factor $e^{i\alpha}$. Further simplifications can be achieved for the axisymmetric problem, see Appendix B5. It is worth noting that the reduced matrices $\hat{\mathbf{M}}_m$ form a block-diagonal matrix

$$ \hat{\mathbf{M}} = \begin{pmatrix} \dots & \dots & \dots & \dots & \dots \\ \dots & \hat{\mathbf{M}}_{-1} & 0 & 0 & \dots \\ \dots & 0 & \hat{\mathbf{M}}_0 & 0 & \dots \\ \dots & 0 & 0 & \hat{\mathbf{M}}_1 & \dots \\ \dots & \dots & \dots & \dots & \dots \end{pmatrix}, \quad (29) $$

which can be seen as a transformation of the original matrix $\mathbf{M}$ by re-ordering its columns and rows. Since these two matrices have the same eigenvalues, it is much faster to search for the eigenvalues of the matrix $\hat{\mathbf{M}}$, which are obtained by combining the eigenvalues of its diagonal blocks $\hat{\mathbf{M}}_m$.

To illustrate the fast convergence of the numerical method, we choose the prolate spheroid with semi-axes $a = 0.5$ and $b = 1$ and compute three eigenvalues $\mu_{00}(n_{max})$, $\mu_{01}(n_{max})$, and $\mu_{02}(n_{max})$ of the truncated matrix $\mathbf{M}$ of size $(n_{max}+1)^2 \times (n_{max}+1)^2$ as functions of the truncation order $n_{max}$. To estimate the truncation error, we subtract from $\mu_{0n}(n_{max})$ its value at $n_{max} = 20$, considered as a proxy of the limiting value $\mu_{0n}(\infty)$. Figure 2 shows how fast the error decreases with $n_{max}$. We checked the rapid convergence for other values of $m$ and $n$, as well as for various aspect ratios (not shown). In particular, we found that larger truncation orders may be needed to achieve high accuracy when $a/b$ gets smaller.

<figure><img src="image_4.png" /><figcaption>FIG. 2. Illustration for the convergence of the numerical method in the case of the exterior of the prolate spheroid with semi-axes $a = 0.5$ and $b = 1$. Three eigenvalues $\mu_{00}(n_{max})$, $\mu_{01}(n_{max})$, and $\mu_{02}(n_{max})$ of the truncated matrix $\mathbf{M}$ of size $(n_{max}+1)^2 \times (n_{max}+1)^2$ as functions of the truncation order $n_{max}$ are shown by symbols. Their values at $n_{max} = 20$, considered here as a benchmark, are subtracted to estimate the error of truncation.</figcaption></figure>

Moreover, when $a/b$ is close to 1, the matrix elements $F_{n,n'}^{m}(\cosh \alpha_0)$ vanish according to Eq. (26) that requires a suitable rescaling and further improvements of the numerical procedure described in Appendix B. Since the case of a sphere is fully explicit, we do not study geometric settings of almost spherical domains.

A similar error analysis was performed for other problems (interior of a prolate spheroid, oblate spheroids, see below). In the following, we generally use $n_{max} = 10$ that is sufficient to produce accurate numerical results for all the considered settings.

### Examples of eigenfunctions

Figure 3 illustrates ten eigenfunctions $v_{mn}$ of the Dirichlet-to-Neumann operator $\mathcal{M}$ for the exterior of the prolate spheroid with semi-axes $a = 0.5$ and $b = 1$. The ground eigenfunction $v_{00}$ is not constant (see below), even though its minor changes are difficult to see due to the chosen colorbar, for which color changes in the same range of values from $-1$ to $1$ for all shown eigenfunctions. The geometric structure of the remaining shown eigenfunctions resembles that of the spherical harmonics $Y_{mn}(\theta, \phi)$. Note that the eigenfunctions $v_{mn}$ and $v_{(-m)n}$ correspond to the same eigenvalue $\mu_{mn}$ and differ only by the factor $e^{\pm im\phi}$; for this reason, the eigenfunctions $v_{(-m)n}$ are not shown.

Figure 4 presents the behavior of $v_{mn}(\theta, \phi)$ as a function of $\theta$ for $\phi = 0$ (i.e., its projection onto the $xz$ plane) for prolate spheroidal surfaces of variables minor semi-axis $a$ (with fixed $b$). As $a \to b$, the surface becomes spherical, and the eigenfunctions $v_{mn}$ coincide with normalized spherical harmonics $Y_{mn}$. In turn, as $a$ decreases,
<figure><img src="image_5.png" /><figcaption>FIG. 3. Several eigenfunctions $v_{mn}$ of the Dirichlet-to-Neumann operator $\mathcal{M}$ for the exterior of the prolate spheroid with semi-axes $a = 0.5$ and $b = 1$. The associated eigenvalues are shown on the top. The eigenfunctions with even $m + n$ are symmetric with respect to the horizontal plane $z = 0$, whereas the eigenfunctions with odd $m + n$ are antisymmetric, in agreement with Eq. (32). The truncation order is $n_{max} = 10$.</figcaption></figure>

<table><caption>TABLE I. Several eigenvalues $\mu_{mn}$ of the Dirichlet-to-Neumann operator $\mathcal{M}$ for the exterior of the prolate spheroid with semi-axes $a$ and $b = 1$. The matrices $\hat{M}_m$ are truncated to the size $(n_{max} + 1 - m) \times (n_{max} + 1 - m)$, with $n_{max} = 10$. Further increase of $n_{max}$ did not change the shown eigenvalues.</caption><thead><tr><th>a/b</th><th>μ<sub>00</sub></th><th>μ<sub>01</sub></th><th>μ<sub>02</sub></th><th>μ<sub>11</sub></th><th>μ<sub>12</sub></th><th>μ<sub>22</sub></th></tr></thead><tbody><tr><td>0.1</td><td>3.960</td><td>6.787</td><td>8.740</td><td>11.041</td><td>12.764</td><td>20.748</td></tr><tr><td>0.2</td><td>2.558</td><td>4.820</td><td>6.523</td><td>6.035</td><td>7.649</td><td>10.797</td></tr><tr><td>0.3</td><td>2.019</td><td>4.011</td><td>5.577</td><td>4.379</td><td>5.924</td><td>7.512</td></tr><tr><td>0.4</td><td>1.717</td><td>3.514</td><td>4.964</td><td>3.554</td><td>5.034</td><td>5.889</td></tr><tr><td>0.5</td><td>1.516</td><td>3.145</td><td>4.489</td><td>3.057</td><td>4.467</td><td>4.926</td></tr><tr><td>0.6</td><td>1.367</td><td>2.844</td><td>4.092</td><td>2.721</td><td>4.054</td><td>4.287</td></tr><tr><td>0.7</td><td>1.250</td><td>2.588</td><td>3.754</td><td>2.475</td><td>3.726</td><td>3.832</td></tr><tr><td>0.8</td><td>1.153</td><td>2.365</td><td>3.464</td><td>2.285</td><td>3.451</td><td>3.490</td></tr><tr><td>0.9</td><td>1.071</td><td>2.171</td><td>3.215</td><td>2.130</td><td>3.212</td><td>3.220</td></tr><tr><td>1.0</td><td>1</td><td>2</td><td>3</td><td>2</td><td>3</td><td>3</td></tr></tbody></table>

the eigenfunctions $v_{mn}$ deviate further and further from spherical harmonics. Note that $v_{00}$ and $v_{01}$ are axisymmetric (independent of $\phi$), so that their structure, shown for $\phi = 0$, remains the same for any $\phi$. In turn, the structure of $v_{11}$ is affected by the factor $e^{i\phi}$; in particular, it is complex-valued, and its real part becomes negative for $\phi$ between $\pi/2$ and $3\pi/2$. Table I summarizes the first eigenvalues $\mu_{mn}$ for various prolate spheroids.

<figure><img src="image_6.png" /><figcaption>FIG. 4. The $xz$ projection of the eigenfunctions $v_{mn}(\theta, \phi)$ of the Dirichlet-to-Neumann operator $\mathcal{M}$ for the exterior of the prolate spheroid with semi-axis $a$ and $b = 1$, with $mn = 00$ (top), $mn = 01$ (middle), and $mn = 11$ (bottom). Thin line shows the normalized spherical harmonics $Y_{mn}(\theta, 0)$ that corresponds to $a/b = 1$. The truncation order is $n_{max} = 10$.</figcaption></figure>

## Asymptotic behavior for elongated spheroids

In the limit $a \to 0$, prolate spheroids get thinner and thinner, approaching a needle of length $2b$. In this limit, one has $\alpha_0 \to 0$, and the coefficients $c_{mn}$ diverge in the leading order as

$$c_{mn} \approx \begin{cases} \frac{1}{a \ln(b/a)} & (m = 0), \\ \frac{|m|}{a} & (m \neq 0). \end{cases} \qquad (30)$$
<figure><img src="image_7.png" /><figcaption>FIG. 5. Several eigenvalues $\mu_{mn}$ of the Dirichlet-to-Neumann operator $\mathcal{M}$ for the exterior of the prolate spheroid as functions of its minor semi-axis $a$, with $b = 1$. At $a = b$, one retrieves the eigenvalues $(n+1)/b$ for the exterior of a sphere of radius $b$. Thin lines from bottom to top present $1/(a \ln(b/a))$, $1/a$ and $2/a$ that capture the asymptotic behavior of $\mu_{mn}$ with $m = 0, 1, 2$, respectively. The eigenvalues are obtained by diagonalizing $\mathbf{M}_m$ truncated to the size $(n_{\max} + 1 - m) \times (n_{\max} + 1 - m)$, with $n_{\max} = 10$.</figcaption></figure>

In turn, the coefficients $F_{n,n'}^m(1)$ are finite. As a consequence, the eigenvalues $\mu_{mn}$ of the matrix $\mathbf{M}_m$ diverge as $a \to 0$ as

$$ \mu_{mn} \approx \begin{cases} \frac{q_n}{a \ln(b/a)} (1 + O(1/\ln(b/a))) & (m=0), \\ \frac{|m|}{a} & (m \neq 0), \end{cases} \quad (31) $$

with some prefactors $q_n$. Figure 5 illustrates the behavior of the first eigenvalues $\mu_{mn}$ as functions of the minor semi-axis $a$ (with fixed $b=1$). At $a=1$, one retrieves the eigenvalues $(n+1)/b$ for the exterior of a sphere of radius $b$. In turn, the eigenvalues $\mu_{mn}$ diverge as $a \to 0$ according to Eq. (31). Note that the asymptotic behavior (31) with the numerical prefactor $q_0 \approx 1$ is quite accurate for $\mu_{00}$. In turn, logarithmic corrections to the leading order are more significant for $\mu_{0n}$ with $n > 0$.

### Reflection symmetry

Since the domain $\Omega$ is symmetric under reflection with respect to the horizontal plane $z=0$, the Steklov eigenfunctions inherit this symmetry. In fact, if $V(x, y, z)$ is a Steklov eigenfunction, corresponding to an eigenvalue $\mu$, then $V(x, y, -z)$ is also an eigenfunction corresponding to the same $\mu$. Moreover, their linear combinations $V_{\pm}(x, y, z) = V(x, y, z) \pm V(x, y, -z)$ are also eigenfunctions, if they are not zero. There are thus three options: (i) $V(x, y, z)$ is symmetric: $V(x, y, z) = V(x, y, -z)$; (ii) $V(x, y, z)$ is antisymmetric: $V(x, y, z) = -V(x, y, -z)$; or (iii) $V(x, y, z)$ is neither symmetric, nor antisymmetric, in which case $V_+(x, y, z)$ and $V_-(x, y, z)$ are symmetric and antisymmetric, respectively.

The structure of the matrices $\mathbf{M}_m$ implies that the Steklov eigenfunctions $V_{mn}(\alpha, \theta, \phi)$, determined by Eq. (24b), are symmetric (resp., antisymmetric) under reflection with respect to the horizontal plane $z=0$ when $m+n$ is even (resp., odd). Indeed, the relation (25) causes that the matrix elements $\mathbf{M}_{mn,mn'}$ are zero when $n+n'$ is odd. As a consequence, the elements $[\mathbf{V}_{mn}]_{mn'}$ of its eigenvector $\mathbf{V}_{mn}$ are zero when $n+n'$ is odd, so that

$$ v_{mn}(\pi - \theta, \phi) = \sum_{n'=0}^{\infty} [\mathbf{V}_{mn}]_{mn'} (-1)^{m+n'} Y_{mn'}(\theta, \phi), $$

where we used $Y_{mn'}(\pi - \theta, \phi) = (-1)^{m+n'} Y_{mn'}(\theta, \phi)$ according to Eq. (25). As the terms with odd $n+n'$ vanish, one concludes that

$$ v_{mn}(\alpha, \pi - \theta, \phi) = (-1)^{m+n} v_{mn}(\theta, \phi). \quad (32) $$

This reflection symmetry is clearly illustrated in Fig. 3. The same symmetry is preserved for the Steklov eigenfunctions:

$$ V_{mn}(\alpha, \pi - \theta, \phi) = (-1)^{m+n} V_{mn}(\alpha, \theta, \phi). \quad (33) $$

### Half-spheroid in the half-space

The above symmetries provide a complementary insight onto the Steklov problem in the upper half-space. Let

$$ \begin{aligned} \Omega_+ &= \left\{ (x, y, z) \in \mathbb{R}^3 : \frac{x^2}{a^2} + \frac{y^2}{a^2} + \frac{z^2}{b^2} > 1, z > 0 \right\} \\ &= \left\{ (\alpha, \theta, \phi) : \alpha > \alpha_0, 0 \le \theta < \pi/2, 0 \le \phi < 2\pi \right\} \end{aligned} $$

be the exterior of a prolate spheroid in the upper half-space (the second line highlights that the upper half-space corresponds to $0 \le \theta < \pi/2$ in the prolate spheroidal coordinates). The boundary of this domain is the union of the upper half-spheroidal surface,

$$ \partial\Omega_+ = \{ (\alpha, \theta, \phi) : \alpha = \alpha_0, 0 \le \theta < \pi/2, 0 \le \phi < 2\pi \}, $$

and the remaining horizontal plane at $z=0$ with a circular hole of radius $a$:

$$ \partial\Omega_0 = \{ (\alpha, \theta, \phi) : \alpha > \alpha_0, \theta = \pi/2, 0 \le \phi < 2\pi \}. $$

According to Eq. (33), antisymmetric Steklov eigenfunctions $V_{mn}$ (with odd $m+n$) vanish at $\theta = \pi/2$ that corresponds to the Dirichlet boundary condition on $\partial\Omega_0$. As a consequence, they solve the mixed Steklov-Dirichlet exterior problem:

$$ \Delta V_{mn} = 0 \quad \text{in } \Omega_+, \quad (34a) $$

$$ \partial_n V_{mn} = \mu_{mn} V_{mn} \quad \text{on } \partial\Omega_+, \quad (34b) $$

$$ V_{mn} = 0 \quad \text{on } \partial\Omega_0 \quad (34c) $$
(here the notation $\partial_n$ for the normal derivative should not be confused with the index $n$). Equivalently, one can speak of eigenvalues $\mu_{mn}$ and eigenfunctions $v_{mn}$ of the Dirichlet-to-Neumann operator $M^D$ that maps a given function $f$ on the half-spheroidal boundary $\partial\Omega_+$ onto another function $g = M^D f = (\partial_n u)|_{\partial\Omega_+}$, where $u$ satisfies

$$ \Delta u = 0 \quad \text{in } \Omega_+, \quad u|_{\partial\Omega_+} = f, \quad u|_{\partial\Omega_0} = 0. \quad (35) $$

For instance, $v_{01}$, $v_{03}$, $v_{12}$ and $v_{23}$ shown in Fig. 3 are examples of eigenfunctions of $M^D$.

In turn, the symmetric Steklov eigenfunctions $V_{mn}$ (with even $m+n$) satisfy

$$ (\partial_n V_{mn})|_{\partial\Omega_0} = \left( \frac{1}{h_\theta} \partial_\theta V_{mn} \right)|_{\theta=\pi/2} = 0 $$

that corresponds to the Neumann boundary condition on $\partial\Omega_0$. In other words, they solve the mixed Steklov-Neumann exterior problem:

$$ \Delta V_{mn} = 0 \quad \text{in } \Omega_+, \quad (36a) $$

$$ \partial_n V_{mn} = \mu_{mn} V_{mn} \quad \text{on } \partial\Omega_+, \quad (36b) $$

$$ \partial_n V_{mn} = 0 \quad \text{on } \partial\Omega_0. \quad (36c) $$

Equivalently, one can speak of eigenvalues $\mu_{mn}$ and eigenfunctions $v_{mn}$ of the Dirichlet-to-Neumann operator $M^N$ that maps a given function $f$ on the half-spheroidal boundary $\partial\Omega_+$ onto another function $g = M^N f = (\partial_n u)|_{\partial\Omega_+}$, where $u$ satisfies

$$ \Delta u = 0 \quad \text{in } \Omega_+, \quad u|_{\partial\Omega_+} = f, \quad (\partial_n u)|_{\partial\Omega_0} = 0. \quad (37) $$

For instance, $v_{00}$, $v_{02}$, $v_{11}$, $v_{22}$, $v_{13}$ and $v_{33}$ shown in Fig. 3 are examples of eigenfunctions of $M^N$.

### III. OBLATE SPHEROIDS

For the exterior of an oblate spheroid with semi-axes $a \le b$,

$$ \Omega = \left\{ (x, y, z) \in \mathbb{R}^3 : \frac{x^2}{b^2} + \frac{y^2}{b^2} + \frac{z^2}{a^2} > 1 \right\}, \quad (38) $$

the computation is very similar so that we only sketch the main steps and formulas. In the oblate spheroidal coordinates $(\alpha, \theta, \phi)$,

$$ \begin{pmatrix} x \\ y \\ z \end{pmatrix} = a_E \begin{pmatrix} \cosh \alpha \cos \theta \cos \phi \\ \cosh \alpha \cos \theta \sin \phi \\ \sinh \alpha \sin \theta \end{pmatrix} \quad \left\{ \begin{array}{l} 0 \le \alpha < \infty \\ -\frac{\pi}{2} \le \theta \le \frac{\pi}{2} \\ 0 \le \phi < 2\pi \end{array} \right\}, $$

with $a_E = \sqrt{b^2 - a^2}$ (Fig. 1b), the scale factors determining the surface and volume elements, are [43]

$$ h_\alpha = h_\theta = a_E \sqrt{\sinh^2 \alpha + \sin^2 \theta}, $$

$$ h_\phi = a_E \cosh \alpha \cos \theta. $$

In these coordinates, the domain $\Omega$ is still characterized by $\alpha > \alpha_0 = \tanh^{-1}(a/b)$, and the action of the Laplace operator reads

$$ \Delta u = \frac{1}{a_E^2(\sinh^2 \alpha + \sin^2 \theta)} \left[ \frac{1}{\cosh \alpha} \frac{\partial}{\partial \alpha} \left( \cosh \alpha \frac{\partial u}{\partial \alpha} \right) + \frac{1}{\cos \theta} \frac{\partial}{\partial \theta} \left( \cos \theta \frac{\partial u}{\partial \theta} \right) \right] + \frac{1}{a_E^2 \cosh^2 \alpha \cos^2 \theta} \frac{\partial^2 u}{\partial \phi^2}. \quad (39) $$

A general solution of the Laplace equation reads

$$ u(\alpha, \theta, \phi) = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} A_{mn} Q_n^m(i \sinh \alpha) \bar{Y}_{mn}(\theta, \phi), \quad (40) $$

where $\bar{Y}_{mn}(\theta, \phi) = a_{mn} P_n^m(\sin \theta) e^{im\phi}$, and the action of the Dirichlet-to-Neumann operator is

$$ \begin{aligned} \mathcal{M}u|_{\partial\Omega} &= \left.\frac{\partial u}{\partial n}\right|_{\partial\Omega} = -\frac{1}{h_\alpha} \left.\partial_\alpha u\right|_{\alpha=\alpha_0} \quad (41) \\ &= \sum_{n=0}^{\infty} \sum_{m=-n}^{n} A_{mn} \bar{Y}_{mn}(\theta, \phi) \frac{\cosh \alpha_0 Q_n^{m'}(i \sinh \alpha_0)}{i a_E \sqrt{\cosh^2 \alpha_0 - \cos^2 \theta}}. \end{aligned} $$

Multiplying this relation by $\bar{Y}_{m'n'}^*(\theta, \phi) \cos \theta$ and integrating over $\theta$ and $\phi$, one gets a matrix representation of the operator $\mathcal{M}$ on the orthonormal basis of $\bar{Y}_{mn}$:

$$ \mathbf{M}_{mn,m'n'} = 2\pi\delta_{m,m'}a_{mn}a_{m'n'}c_{m'n'}\bar{F}_{n,n'}^{m}(\sinh \alpha_0), \quad (42) $$

where

$$ c_{mn} = \frac{\cosh \alpha_0 Q_n^{m'}(i \sinh \alpha_0)}{i a_E Q_n^m(i \sinh \alpha_0)} \quad (43) $$

and

$$ \bar{F}_{n,n'}(z) = \int_{-1}^{1} dx \frac{P_{n}^{m}(x) P_{n'}^{m'}(x)}{\sqrt{z^2 + x^2}}. \quad (44) $$

According to Eq. (15), one also has

$$ \bar{F}_{n,n'}^{m}(z) = i F_{n,n'}^{m}(iz). \quad (45) $$

Since the structure of the matrix $\mathcal{M}$ is the same as for prolate spheroids, many properties of eigenfunctions of the Dirichlet-to-Neumann operator remain unchanged; in particular, they can be classified according to their dependence on $\phi$ via $e^{im\phi}$; we employ the double index $mn$ in the following. Using again the decomposition (22), one can diagonalize separately the matrices $\mathcal{M}_m$ to access the eigenvalues $\mu_{mn}$ and to construct the eigenfunctions

$$ v_{mn}(\theta, \phi) = \sum_{n'=0}^{\infty} [\mathbf{V}_{mn}]_{mn'} \bar{Y}_{mn'}(\theta, \phi) \quad (46) $$

and the Steklov eigenfunctions

$$ V_{mn}(\alpha, \theta, \phi) = \sum_{n'=0}^{\infty} \frac{Q_{n'}^m(i \sinh \alpha)}{Q_n^m(i \sinh \alpha_0)} [\mathbf{V}_{mn}]_{mn'} \bar{Y}_{mn'}(\theta, \phi). \quad (47) $$
<figure><img src="image_8.png" /><figcaption>FIG. 6. Several eigenfunctions $v_{mn}$ of the Dirichlet-to-Neumann operator $\mathcal{M}$ for the exterior of the oblate spheroid with semi-axes $a = 0.5$ and $b = 1$. The associated eigenvalues are shown on the top. The eigenfunctions with even $m+n$ are symmetric (with respect to the horizontal plane $z = 0$), whereas the eigenfunctions with odd $m+n$ are antisymmetric, in agreement with Eq. (32). The truncation order is $n_{max} = 10$.</figcaption></figure>

As previously, the Steklov eigenfunctions $V_{mn}$ are symmetric for even $m+n$ and antisymmetric for odd $m+n$ under reflection with respect to the horizontal plane $z = 0$. In particular, these eigenfunctions solve the mixed Steklov-Neumann and Steklov-Dirichlet exterior problems for the exterior of an oblate spheroid in the upper half-space.

Figure 6 illustrates several eigenfunctions $v_{mn}$ of the Dirichlet-to-Neumann operator $\mathcal{M}$ for the exterior of the oblate spheroid with semi-axes $a = 0.5$ and $b = 1$. Expectedly, these eigenfunctions resemble spherical harmonics $\bar{Y}_{mn}$ but exhibit some differences. In particular, the ground eigenfunction $v_{00}$ is not constant, though its variations are small.

Figure 7 illustrates the behavior of three eigenvalues $\mu_{0n}$ as functions of the minor semi-axis $a$ (with fixed $b = 1$). At $a = 1$, one retrieves the eigenvalues $(n+1)/b$ for the exterior of a sphere of radius $b$. As $a$ decreases, oblate spheroids become thinner and thinner, approaching a disk of radius $b$. In contrast to the case of prolate spheroids, the eigenvalues are finite in this limit. Curiously, each eigenvalue does not change monotonously with $a$. Table II summarizes several eigenvalues $\mu_{mn}$ for oblate spheroids with the minor semi-axis $a$ ranging from 0 to 1 (with $b = 1$).

## Limit of a disk

Let us focus on the limit of the disk when $a = 0$ and thus $\alpha_0 = 0$. One can introduce the radial coordinate is

<figure><img src="image_9.png" /><figcaption>FIG. 7. Three eigenvalues $\mu_{0n}$ of the Dirichlet-to-Neumann operator for the exterior of an oblate spheroid as functions of its minor semi-axis $a$, with $b = 1$. At $a = b$, one retrieves the eigenvalues $(n+1)/b$ for the exterior of a sphere of radius $b$. As $a \to 0$, one gets the eigenvalues for the exterior of a disk of radius $b$. The axisymmetric matrix $\mathbf{M}_0$ from Eq. (B44) was truncated to the size $(n_{max}+1) \times (n_{max}+1)$, with $n_{max} = 20$.</figcaption></figure>

<table><caption>TABLE II. Several eigenvalues $\mu_{mn}$ of the Dirichlet-to-Neumann operator $\mathcal{M}$ for the exterior of the oblate spheroid with semi-axes $a$ and $b = 1$. The matrices $\hat{M}_m$ are truncated to the size $(n_{max}+1-m) \times (n_{max}+1-m)$, with the truncation order $n_{max} = 10$. Further increase of $n_{max}$ or decrease of $\delta$ do not change the shown eigenvalues.</caption><thead><tr><th>a/b</th><th>μ<sub>00</sub></th><th>μ<sub>01</sub></th><th>μ<sub>02</sub></th><th>μ<sub>11</sub></th><th>μ<sub>12</sub></th><th>μ<sub>22</sub></th></tr></thead><tbody><tr><td>0.0</td><td>1.158</td><td>2.006</td><td>4.317</td><td>2.755</td><td>3.453</td><td>4.121</td></tr><tr><td>0.1</td><td>1.204</td><td>2.057</td><td>4.314</td><td>2.811</td><td>3.512</td><td>4.197</td></tr><tr><td>0.2</td><td>1.220</td><td>2.094</td><td>4.206</td><td>2.796</td><td>3.539</td><td>4.166</td></tr><tr><td>0.3</td><td>1.217</td><td>2.117</td><td>4.040</td><td>2.732</td><td>3.536</td><td>4.058</td></tr><tr><td>0.4</td><td>1.200</td><td>2.129</td><td>3.850</td><td>2.634</td><td>3.506</td><td>3.900</td></tr><tr><td>0.5</td><td>1.173</td><td>2.128</td><td>3.666</td><td>2.520</td><td>3.453</td><td>3.721</td></tr><tr><td>0.6</td><td>1.141</td><td>2.118</td><td>3.498</td><td>2.401</td><td>3.381</td><td>3.543</td></tr><tr><td>0.7</td><td>1.106</td><td>2.098</td><td>3.350</td><td>2.286</td><td>3.295</td><td>3.379</td></tr><tr><td>0.8</td><td>1.070</td><td>2.071</td><td>3.220</td><td>2.180</td><td>3.200</td><td>3.234</td></tr><tr><td>0.9</td><td>1.034</td><td>2.038</td><td>3.105</td><td>2.085</td><td>3.101</td><td>3.108</td></tr><tr><td>1.0</td><td>1</td><td>2</td><td>3</td><td>2</td><td>3</td><td>3</td></tr></tbody></table>

$r = \sqrt{x^2 + y^2} = b|\cos\theta|$, with $0 < \theta < \pi/2$ corresponding to the upper side of the disk, and $-\pi/2 < \theta < 0$ corresponding to its lower side. Formally, one can set $r = b \cos\theta$ and associate positive $r$ with the upper side of the disk and negative $r$ with its lower side.

As $\alpha_0 \to 0$, the functions $\hat{F}_{n,n'}^m(z)$ from Eq. (44) logarithmically diverge in the limit $z = \sinh\alpha_0 \to 0$. For any $\alpha_0 > 0$, one can still use the truncated matrix $\mathbf{M}$ (or $\mathbf{M}_m$) to approximate the eigenvalues and eigenfunctions but larger and larger truncation orders are needed as $\alpha_0$ decreases. For this reason, it is convenient to use an alternative matrix representation, which corresponds to the operator $\mathcal{M}^{-1}$ and remains valid even in the limit $\alpha_0 = 0$. This representation is described in Appendix
<figure><img src="image_10.png" /><figcaption>FIG. 8. Several eigenfunctions $v_{mn}$ of the Dirichlet-to-Neumann operator $\mathcal{M}$ for the exterior of the disk of radius $b = 1$. The associated eigenvalues are shown on the top. Normalization by the maximum of $|v_{mn}|$ was employed for a better visualization. The symmetric (resp., antisymmetric) shape of the eigenfunction $v_{mn}$ for even (resp., odd) $m+n$ allows one to reconstruct its structure on the bottom side of the disk, which is not visible. The truncation order is $n_{max} = 10$.</figcaption></figure>

## B 2.

Figure 8 illustrates several eigenfunctions $v_{mn}$ of the Dirichlet-to-Neumann operator for the exterior of the disk of radius $b = 1$. One can see that the ground eigenfunction $v_{00}$ shows more significant variations as compared to that shown in Fig. 6. We stress that the eigenfunctions are also present on the bottom side of the disk, which is not visible. Their structure on this hidden side can be easily reconstructed from their symmetry with respect to the horizontal plane: $v_{mn}$ is symmetric (resp., antisymmetric) for even (resp., odd) $m+n$. For instance, the eigenfunction $v_{01}$, which is positive on the upper side of the disk, takes negative values on the bottom side.

### Axisymmetric setting

Let us briefly discuss the action of the Dirichlet-to-Neumann operator onto rotationally invariant functions that do not depend on the angle $\phi$. If a function $f(\theta)$ is decomposed on the normalized Legendre polynomials $\psi_n(\theta) = \sqrt{n+1/2} P_n(\sin \theta)$,

$$f(\theta) = \sum_{n=0}^{\infty} f_n \psi_n(\theta), \qquad (48)$$

the action of $\mathcal{M}$ reads

$$\mathcal{M}f = \sum_{n=0}^{\infty} f_n \psi_n(\theta) \frac{c_{0n}}{b\sqrt{1 - \cos^2 \theta}}. \qquad (49)$$

Setting $b = 1$ for simplicity and using $r = \cos \theta$, we get, for instance,

$$\mathcal{M}1 = \frac{c_{00}}{\sqrt{1 - r^2}}, \qquad (50a)$$

$$\mathcal{M}\sqrt{1 - r^2} = c_{01}, \qquad (50b)$$

<figure><img src="image_11.png" /><figcaption>FIG. 9. Four eigenfunctions $v_{0n}(\theta)$ of the Dirichlet-to-Neumann operator for the exterior of a disk of radius $b = 1$, as a function of $r = b \cos \theta$. The associated eigenvalues are 1.158, 2.006, 4.317, and 5.125. The axisymmetric matrix $\mathcal{M}_0$ is truncated to $(n_{max} + 1) \times (n_{max} + 1)$, with $n_{max} = 20$.</figcaption></figure>

where

$$c_{0n} = \frac{-iQ'_{n}(0)}{Q_{n}(0)} = 2\left(\frac{\Gamma(n/2 + 1)}{\Gamma(n/2 + 1/2)}\right)^{2} \qquad (51)$$

(see Appendix B5 for details). For instance, $c_{00} = 2/\pi$, $c_{01} = \pi/2$, $c_{02} = 8/\pi$, $c_{03} = 9\pi/8$, etc. The first relation (50a) reproduces the classical Weber's solution for the electric current density onto a conducting disk [44], see below.

Alternatively, setting $x = \sin \theta$, one has

$$\mathcal{M}P_n(x) = c_{0n} \frac{P_n(x)}{|x|} \qquad (-1 < x < 1). \qquad (52)$$

One sees that $\mathcal{M}P_n(x)$ keeps the parity of $P_n(x)$, i.e., it is symmetric for even $n$ and antisymmetric for odd $n$. For instance, one has

$$\mathcal{M}1 = \frac{2}{\pi} \frac{1}{|x|},$$

$$\mathcal{M}x = \frac{\pi}{2} \frac{x}{|x|},$$

$$\mathcal{M}x^2 = \frac{2}{\pi} \frac{4x^2 - 1}{|x|}.$$

### Disk in the half-space

The general representation (40) allows one to solve boundary value problems (35, 37) in the upper half-space, by imposing Dirichlet or Neumann boundary condition on $\partial \Omega_0$. In the case of a disk of radius $b$, the first problem actually concerns the Dirichlet boundary condition on the horizontal plane, for a class of functions $f(r)$ that
are strictly zero for $r \ge b$. Its solution can be written in the integral form:

$$u(x_0, y_0, z_0) = \int_{x^2+y^2<b^2} dx dy f(x, y) \times \frac{z_0}{2\pi(z_0^2+(x-x_0)^2+(y-y_0)^2)^{3/2}}, \quad (53)$$

where the factor in front of $f$ is the harmonic measure density [45]. Alternatively, one can use the representation (40) of a harmonic function in oblate spheroidal coordinates, in which the coefficients $A_{mn}$ are obtained by setting $\alpha = 0$, multiplying by $\bar{Y}_{mn}^* \cos \theta$ and integrating over $\theta$ and $\phi$:

$$A_{mn} = \frac{1}{Q_n^m(0)} \int_{-\pi/2}^{\pi/2} d\theta \int_0^{2\pi} d\phi f(\theta, \phi) \bar{Y}_{mn}^*(\theta, \phi) \cos \theta. \quad (54)$$

For instance, if $f(\theta, \phi) = \text{sign}(\theta)$, one gets

$$A_{mn} = \frac{2\pi\delta_{m,0}a_{0n}}{Q_n^m(0)} \int_{-1}^{1} dx \sinh(x) P_n(x) \\ = \frac{2\pi\delta_{m,0}a_{0n}}{Q_n^m(0)} \frac{P_{n-1}(0) - P_{n+1}(0)}{n+1/2} \quad (55)$$

for odd $n$, where $P_{2n}(0) = \frac{(-1)^n(2n)!}{2^{2n}(n!)^2}$. As a consequence, one has

$$u = \sum_{n=0}^{\infty} \frac{Q_{2n+1}(i \sinh \alpha)}{Q_{2n+1}(0)} P_{2n+1}(\sin \theta) (P_{2n}(0) - P_{2n+2}(0)). \quad (56)$$

This is also the solution of the Laplace equation in the upper half-space with the boundary condition $u|_{z=0} = \Theta(b-r)$, where $\Theta(z)$ is the Heaviside step function: $\Theta(z) = 1$ for $z > 0$ and $0$ otherwise. This is the splitting probability for the disk of radius $b$, i.e., the probability that a particle started from a point $(\alpha, \theta, \phi)$ hits the disk before hitting its complement (the horizontal plane without this disk) [45–47].

In the same vein, the representation (40) with even $m+n$ allows one to solve the mixed Dirichlet-Neumann boundary value problem (37) in the upper half-space, with coefficients $A_{mn}$ given by Eq. (54). For instance, if $f=1$, one gets

$$A_{mn} = \frac{2\sqrt{\pi}\delta_{m,0}\delta_{n,0}}{Q_n(0)}, \quad (57)$$

so that

$$u = \frac{Q_0(i \sinh \alpha)}{Q_0(0)}, \quad (58)$$

where $\alpha = \cosh^{-1}((r_+ + r_-)/(2b))$, with $r_\pm = \sqrt{(r \pm b)^2 + z^2}$. This expression can be written more explicitly as

$$u(r, z) = \frac{1}{\pi} \cos^{-1} \left( 1 - \frac{8b^2}{(r_+ + r_-)^2} \right). \quad (59)$$

It is worth noting that the mixed Dirichlet-Neumann boundary value problem for rotationally invariant functions $f(r)$,

$$\Delta u = 0 \quad (z > 0), \quad \begin{cases} u(r, 0) = f(r) & (0 \le r < b), \\ (\partial_z u(r, z))|_{z=0} = 0 & (r > b), \end{cases} \quad (60)$$

was solved by Beltrami via integral representations (see [44], page 66):

$$u(r, z) = \int_{0}^{\infty} d\xi A(\xi) e^{-z\xi} J_{0}(\xi r), \quad (61)$$

where

$$A(\xi) = \int_{0}^{1} ds \cos(\xi s) \frac{d}{ds} \left( \frac{2}{\pi} \int_{0}^{s} dt \frac{tf(tb)}{\sqrt{s^2 - t^2}} \right). \quad (62)$$

For instance, setting $f(r) = 1$, one retrieves the Weber's solution for the electric potential around a conducting disk:

$$u(r, z) = \frac{2}{\pi} \int_{0}^{\infty} \frac{d\xi}{\xi} \sin(\xi b) e^{-\xi z} J_{0}(\xi r), \quad (63)$$

from which

$$-(\partial_n u)|_{\partial \Omega_+} = \frac{2}{\pi \sqrt{b^2 - r^2}} \quad (0 \le r < b), \quad (64)$$

where $J_0(z)$ is the Bessel function of the first kind. One sees that the above expressions (58) and (59) are more explicit than the equivalent integral representation (63). In probabilistic terms, $1 - u(r, z)$ is the escape (or survival) probability for a particle started from a point $(r, z)$, which can be destroyed upon its first arrival onto the absorbing disk. Alternatively, $C_0 u(r, z)$ can be understood as the concentration of particles diffusing from infinity (with a concentration $C_0$) towards an absorbing disk. In this light, $j = -DC_0(\partial_n u)|_{\partial \Omega_+}$ is the diffusive flux density onto the disk. Dividing this expression by the total flux, one gets the harmonic measure density on the absorbing disk for a particle started from infinity:

$$\omega(r) = \frac{DC_0 \frac{2}{\pi \sqrt{b^2 - r^2}}}{DC_0 \frac{2}{\pi} b} = \frac{1}{b \sqrt{b^2 - r^2}}. \quad (65)$$

## IV. TWO APPLICATIONS

As mentioned in Sec. I, the Dirichlet-to-Neumann operator $\mathcal{M}$ and its eigenfunctions have found numerous applications in applied mathematics and engineering. Here we briefly discuss two applications in physics. In Sec. IV A, we start with diffusion-controlled reactions in chemical physics and illustrate the use of $\mathcal{M}$ for representing the Robin boundary condition that is often used
to describe surface reactions on partially reactive targets. In Sec. IV B, we highlight their relation to the statistics of encounters between a diffusing particle and the boundary, which is relevant in statistical physics and in the theory of reflected Brownian motion.

### A. Partially reactive boundaries

Diffusion-controlled reactions play an important role in physics, chemistry and biology [48–58]. In the idealized setting introduced by Smoluchowski [59], a reactant diffuses towards its partner or towards a catalytic surface, and reacts upon their first encounter. In many situations, however, the first encounter is not sufficient, as the reactant and/or its partner may need to overcome an activation energy barrier, to be in an appropriate “active” state, to arrive onto a specific catalytic germ, to pass through a channel/pore, etc. [58]. If any of these conditions is failed, the reactant resumes its bulk diffusion until the next encounter, and so on. Starting from Collins and Kimball [60], such partial reactions are described by imposing the Robin boundary condition, in which the diffusive flux of particles towards the catalytic boundary is postulated to be proportional to their concentration on that boundary. In other words, the concentration C of diffusing reactants in the steady-state regime obeys:

$$ \Delta C = 0 \quad \text{in } \Omega, \quad (66a) $$

$$ -D\partial_n C = \kappa C \quad \text{on } \partial\Omega, \quad (66b) $$

$$ \lim_{|\mathbf{x}| \to \infty} C = C_0, \quad (66c) $$

where D is the diffusion coefficient, $\kappa \ge 0$ is a constant reactivity of the boundary, and $C_0$ is the initial concentration maintained at infinity. Setting $q = \kappa/D$ and $C = C_0(1-u)$, one can express the Robin boundary condition $\partial_n u + qu = q$ on $\partial\Omega$ with the help of the Dirichlet-to-Neumann operator as $M u|_{\partial\Omega} + q u|_{\partial\Omega} = q$. Inverting this relation yields

$$ \begin{aligned} u|_{\partial\Omega} &= (M + q)^{-1}q \\ &= \sum_{k=0}^{\infty} \frac{(v_k, 1)_{L^2(\partial\Omega)}}{\|v_k\|_{L^2(\partial\Omega)}^2} \frac{q}{\mu_k + q} v_k, \end{aligned} \quad (67) $$

where we used the orthogonality and completeness of the basis of the eigenfunctions $\{v_k\}$ of M.

For prolate spheroids, a general solution (7) of the Laplace equation can be written as

$$ u(\alpha, \theta, \phi) = \sum_{m,n} f_{mn} \frac{Q_n^m(\cosh \alpha)}{Q_n^m(\cosh \alpha_0)} Y_{mn}(\theta, \phi), \quad (68) $$

where the coefficients $f_{mn}$ are found from the restriction

of u onto $\partial\Omega$ as:

$$ \begin{aligned} f_{mn} &= \int_0^\pi d\theta \int_0^{2\pi} d\phi \sin\theta Y_{mn}^*(\theta, \phi) u|_{\partial\Omega} \\ &= \sum_{k=0}^\infty \frac{(v_k, 1)_{L^2(\partial\Omega)}}{\|v_k\|_{L^2(\partial\Omega)}^2} \frac{q}{\mu_k + q} \\ &\quad \times \int_0^\pi d\theta \int_0^{2\pi} d\phi \sin\theta Y_{mn}^*(\theta, \phi) v_k(\theta, \phi) \\ &= \sum_{k=0}^\infty \frac{(v_k, 1)_{L^2(\partial\Omega)}}{\|v_k\|_{L^2(\partial\Omega)}^2} \frac{q}{\mu_k + q} [\mathbf{V}_k]_{mn}, \end{aligned} \quad (69) $$

where we employed the representation (17). Using Eqs. (B30, B34) from Appendix B, one can further simplify:

$$ \begin{aligned} f_{mn} &= \frac{\sqrt{4\pi}}{a_E \sinh \alpha_0 Q_0(\cosh \alpha_0)} \sum_{k=0}^{\infty} \frac{[\mathbf{V}_k^*]_{00}[\mathbf{V}_k]_{mn}}{(\mathbf{V}_k^\dagger \mathbf{c} \mathbf{V}_k)} \frac{q}{\mu_k + q} \\ &= \sqrt{4\pi} \sum_{k=0}^{\infty} \frac{[\mathbf{V}_k^*]_{00}[\mathbf{V}_k]_{mn}}{(\mathbf{V}_k^\dagger \mathbf{V}_k)} \frac{q}{\mu_k + q}, \end{aligned} \quad (70) $$

where $\tilde{\mathbf{V}}_k = \mathbf{c}^{\frac{1}{2}} \mathbf{V}_k$, and the diagonal matrix $\mathbf{c}$ is formed by $c_{mn}$. Since $\mu_k$ and $\tilde{\mathbf{V}}_k$ are eigenvalues and eigenvectors of the Hermitian matrix $\tilde{\mathbf{M}} = \mathbf{c}^{\frac{1}{2}} \mathbf{M} \mathbf{c}^{-\frac{1}{2}}$, one can also rewrite the above expression in a matrix form:

$$ f_{mn} = \sqrt{4\pi} q [(\tilde{\mathbf{M}} + q\mathbf{I})^{-1}]_{mn,00}, \quad (71) $$

where $\mathbf{I}$ is the identity matrix. Similar representations were discussed in [17, 40].

Knowing the concentration, one can easily deduce the total diffusive flux onto the boundary:

$$ \begin{aligned} J_q &= \int_{\partial\Omega} (-D\partial_n C) = DC_0 \int_{\partial\Omega} (\partial_n u)|_{\partial\Omega} \\ &= DC_0 q (1, M(M+q)^{-1})_{L^2(\partial\Omega)} \\ &= DC_0 q \sum_{k=0}^{\infty} \frac{|(1, v_k)_{L^2(\partial\Omega)}|^2}{\|v_k\|_{L^2(\partial\Omega)}^2} \frac{\mu_k}{\mu_k + q} \\ &= \frac{4\pi D C_0 q}{\sinh \alpha_0 Q_0^2 (\cosh \alpha_0)} \sum_{k=0}^{\infty} \frac{[\mathbf{V}_k^*]_{00} [\mathbf{V}_k]_{00}}{(\mathbf{V}_k^\dagger \mathbf{c} \mathbf{V}_k)} \frac{1}{\mu_k + q} \\ &= \frac{4\pi D C_0 a_E q}{Q_0 (\cosh \alpha_0)} \sum_{k=0}^{\infty} \frac{[\tilde{\mathbf{V}}_k^*]_{00} [\tilde{\mathbf{V}}_k]_{00}}{(\tilde{\mathbf{V}}_k^\dagger \tilde{\mathbf{V}}_k)} \frac{1}{\mu_k + q} \\ &= \frac{4\pi D C_0 a_E q}{Q_0 (\cosh \alpha_0)} [(\tilde{\mathbf{M}} + q\mathbf{I})^{-1}]_{00,00}. \end{aligned} \quad (72) $$

Denoting by

$$ C = \frac{4\pi a_E}{Q_0(\cosh \alpha_0)} = \frac{8\pi a_E}{\ln\left(\frac{b+a_E}{b-a_E}\right)} \quad (73) $$

the harmonic (Newtonian) capacity of the prolate spheroid, one gets

$$ J_q = J_\infty q [(\tilde{\mathbf{M}} + q\mathbf{I})^{-1}]_{00,00}, \quad (74) $$
where

$$J_{\infty} = DC_{0}C \qquad (75)$$

is the total diffusive flux onto a perfectly reactive prolate spheroid (see also [40]). As a consequence, the effect of partial reactivity is represented by $q[(\tilde{\mathbf{M}}+q\mathbf{I})^{-1}]_{00,00}$. In the limit of a sphere ($a \to b$), this factor reduces to $q/(q+1/b)$ so that one retrieves the Collins-Kimball diffusive flux [60]

$$J_q = 4\pi DC_0 b \frac{qb}{1+qb} \quad (a=b), \qquad (76)$$

where $4\pi DC_0 b$ is the Smoluchowski diffusive flux $J_{\infty}$ onto a perfectly reactive sphere of radius $b$ [59].

For oblate spheroids, a similar computation involves

$$u(\alpha, \theta, \phi) = \sum_{m,n} f_{mn} \frac{Q_n^m(i \sinh \alpha)}{Q_n^m(i \sinh \alpha_0)} \bar{Y}_{mn}(\theta, \phi), \qquad (77)$$

with the same expression (71) for the coefficients $f_{mn}$. The total diffusive flux reads

$$J_q = J_{\infty} q [(\tilde{\mathbf{M}} + q\mathbf{I})^{-1}]_{00,00}, \qquad (78)$$

where

$$J_{\infty} = \frac{4\pi DC_0 a_E}{iQ_0(i \sinh \alpha_0)} = DC_0 C, \qquad (79)$$

and

$$C = \frac{4\pi a_E}{iQ_0(i \sinh \alpha_0)} = \frac{4\pi a_E}{\cos^{-1}(a/b)} \qquad (80)$$

is the harmonic capacity of the oblate spheroid. In the case of a disk, the total flux was thoroughly studied in [61].

Figure 10 illustrates the behavior of the diffusive flux $J_q$, rescaled by $J_{\infty}$, as a function of the reactivity parameter $q$. Expectedly, this ratio changes from 0 at $q = 0$ (no surface reaction) to 1 at $q \to \infty$ (perfect surface reaction). For prolate spheroids, their anisotropy reduces the diffusive flux monotonously (curves are shifted downwards as $a/b$ decreases). As the spheroid gets thinner, its accessibility to Brownian motion is reduced, and this effect is enhanced by decreasing the reactivity $q$ and thus requiring more and more returns to the target to realize a successful reaction. The situation is more subtle for oblate spheroids. Even in the limit $a = 0$, the disk remains accessible to Brownian motion so that the effect of anisotropy onto the diffusive flux is moderate (all shown curves are close to each other). Curiously, the reduction of the diffusive flux by anisotropy is not monotonous for oblate spheroids. For instance, one can notice that the curves for $a/b = 0.01$ and $a/b = 0.2$ cross each other. This non-monotonicity is illustrated on panel (c) which shows the ratio $J_q/J_{\infty}$ as a function of $a/b$. When $q$ is not

<figure><img src="image_12.png" /><figcaption>FIG. 10. The diffusive flux $J_q$ onto a partially reactive spheroidal boundary with semi-axes $a$ and $b$, rescaled by $J_{\infty}$. Panels (a) and (b) present $J_q$ as a function of $qb$ for prolate and oblate spheroids, respectively, whereas panel (c) shows $J_q$ as a function of $a/b$ for oblate spheroids. For an easier comparison between different curves, the ratio $J_q/J_{\infty}$ on panel (c) was rescaled by its minimum (so that all the minimal ratios are equal to 1). The matrix $\tilde{\mathbf{M}}$ in Eqs. (72, 78) was truncated to the size $(n_{\max} + 1)^2 \times (n_{\max} + 1)^2$, with $n_{\max} = 10$.</figcaption></figure>

too large, this ratio exhibits a minimum at some intermediate aspect ratio that depends on the reactivity parameter $q$. From the mathematical point of view, this behavior may be a consequence of the non-monotonous dependence of eigenvalues $\mu_{mn}$ of the Dirichlet-to-Neumann operator on the aspect ratio $a/b$, as illustrated in Fig. 7.
From the physical point of view, it may result from an intricate interplay between the accessibility of the surface controlled by its curvature (i.e., by $a/b$) and the number of failed reaction attempts on it due to partial reactivity (controlled by $qb$). For instance, if the flux $J_q$ was not rescaled by $J_\infty$, there would be no minimum. Further clarifications of this effect require a more detailed analysis of the diffusive dynamics near curved, partially reactive boundaries. We stress, however, that this effect is moderate for oblate spheroids.

More generally, the Dirichlet-to-Neumann operator and its matrix representation allow one to investigate the effect of spatially heterogeneous reactivity when $q$ changes on the boundary [17]. In addition, one can use the Steklov eigenfunctions to compute the Green's function with the Robin boundary condition.

## B. Statistics of boundary encounters

The Steklov eigenfunctions also play the central role in the encounter-based approach to diffusion-controlled reactions [18]. In this context, one usually considers a more general setting, in which the Laplace equation (1) is replaced by the modified Helmholtz equation, $(p - D\Delta)u = 0$, with $p \ge 0$. The eigenfunctions of the related Steklov problem determine the spectral expansion of the Laplace transform of the encounter propagator $P(x, \ell, t|x_0)$, i.e., the joint probability density of the position $X_t$ and the boundary local time $\ell_t$ of a particle that starts from a point $x_0$ and diffuses in a domain $\Omega$ with the reflecting boundary $\partial\Omega$ [18]. Here the boundary local time $\ell_t$ can be understood as the (rescaled) number of encounters between the diffusing particle and the boundary, which is crucial for describing various phenomena occuring on the boundary [20–23]. The encounter propagator allows one to access most common characteristics of diffusion-controlled reactions, including the conventional propagator, the survival probability, the first-passage time probability density, the harmonic measure density, the diffusive flux, etc. For instance, one can investigate the first-crossing time $\tau = \inf\{t > 0 : \ell_t > \ell\}$ of a threshold $\ell$ by the boundary local time $\ell_t$. In other words, the random variable $\tau$ is the first time instance when the number of encounters exceeds a given threshold. The probability density of this random variable,

$$U(\ell, t|x_0)dt = \mathbb{P}_{x_0}\{\tau \in (t, t+dt)\}, \quad (81)$$

was shown to be expressed in terms of the integral of $P(x, \ell, t|x_0)$ over the boundary $\partial\Omega$ (see details in [18, 62]):

$$U(\ell, t|x_0) = D \int_{\partial\Omega} d\mathbf{x} P(\mathbf{x}, \ell, t|x_0). \quad (82)$$

For a bounded domain, any threshold is crossed with probability 1 so that the probability density $U(\ell, t|x_0)$

<figure><img src="image_13.png" /><figcaption>FIG. 11. The crossing probability $P(\ell)$ from Eq. (85) as a function of the threshold $\ell$ for prolate and oblate spheroids with semi-axes $a$ and $b$. The starting point is uniformly distributed on the boundary. The matrix $\tilde{M}$ in Eq. (87) was truncated to the size $(n_{max}+1)^2 \times (n_{max}+1)^2$, with $n_{max} = 10$.</figcaption></figure>

is normalized to 1. However, when a three-dimensional domain is unbounded, the particle can escape to infinity and never return to the boundary $\partial\Omega$. As a consequence, the integral of $U(\ell, t|x_0)$ over $t$ may be smaller than 1; in fact, this integral determines the probability of crossing of a given threshold $\ell$:

$$\mathbb{P}_{x_0}\{\ell_\infty > \ell\} = \int_0^\infty dt U(\ell, t|x_0), \quad (83)$$

which can be seen as the Laplace transform of $U(\ell, t|x_0)$, evaluated at $p = 0$. One can therefore employ the general spectral expansion for the Laplace transform of the encounter propagator $P(x, \ell, t|x_0)$, derived in [18], to determine this probability in terms of the Steklov eigenfunctions for the Laplace equation (i.e., for $p = 0$) that we studied here:

$$\mathbb{P}_{x_0}\{\ell_\infty > \ell\} = \sum_{k=0}^{\infty} \frac{[V_k(x_0)]^*}{\|v_k\|_{L^2(\partial\Omega)}^2} e^{-\mu_k \ell} \int_{\partial\Omega} d\mathbf{x} v_k(\mathbf{x}). \quad (84)$$

At $\ell = 0$, this is simply the probability of hitting the boundary (before escaping to infinity). At large $\ell$, this probability decays exponentially, and the rate of this decay is given by the smallest eigenvalue $\mu_0$. When prolate spheroids get thinner ($a \to 0$), the eigenvalue $\mu_0$ diverges according to Eq. (31), so that the crossing probability $\mathbb{P}_{x_0}\{\ell_\infty > \ell\}$ decays faster with $\ell$. Indeed, it is hard for Brownian motion to access thin prolate spheroids, and it is therefore unlikely to experience many encounters with it (i.e., to get large $\ell_\infty$). In turn, the eigenvalue $\mu_0$ approaches a constant value in the opposite limit of thin oblate spheroids. This is consistent with the fact that even a disk ($a = 0$) remains accessible to Brownian motion in three dimensions.

For illustrative purposes, we consider the starting point to be uniformly distributed over the spheroidal boundary
and evaluate the the associated crossing probability:

$$ 
\begin{aligned}
P(\ell) &= \frac{1}{|\partial\Omega|} \int_{\partial\Omega} d\boldsymbol{x}_0 \mathbb{P}_{\boldsymbol{x}_0}\{\ell_\infty > \ell\} && (85) \\
&= \frac{1}{|\partial\Omega|} \sum_{k=0}^{\infty} e^{-\mu_k \ell} \frac{|(1, v_k)|_{L^2(\partial\Omega)}^2}{\|v_k\|_{L^2(\partial\Omega)}^2}.
\end{aligned}
$$

Comparing this expression with the spectral expansion of the total flux $J_q$ (see the third line in Eq. (72)), one can realize that

$$ \frac{J_q}{D c_0 q |\partial \Omega|} = \int_0^\infty d\ell e^{-q\ell} \rho(\ell), \quad (86) $$

where $\rho(\ell) = -\partial_\ell P(\ell)$ is the probability density of the random variable $\ell_\infty$. As a consequence, we can formally invert the matrix representation (72) to get

$$ P(\ell) = \frac{C}{|\partial\Omega|} [e^{-\ell \tilde{M}} / \tilde{M}]_{00,00} \quad (87) $$

for both prolate and oblate spheroids, with $C$ being their harmonic capacity given by Eqs. (73, 80).

Figure 11 shows the behavior of the crossing probability $P(\ell)$ for prolate and oblate spheroids, as well as for the sphere of radius $b$, for which $P(\ell) = e^{-\ell/b}$. The slowest decay of $P(\ell)$ corresponds to the sphere, which is the most accessible to Brownian motion. As $a/b$ goes to 0, the decay is getting faster and faster for prolate spheroids. In turn, for oblate spheroids, the crossing probability $P(\ell)$ approaches that of a disk. The role of the starting point $\boldsymbol{x}_0$ and other properties of the boundary local time $\ell_\infty$ can be investigated using our results.

V. CONCLUSIONS

In this paper, we investigated the spectral properties of the Dirichlet-to-Neumann operator $\mathcal{M}$ and the related Steklov problem for anisotropic domains in three dimensions. Using prolate and oblate spheroidal coordinates, we derived a matrix representation of $\mathcal{M}$ on the basis of spherical harmonics for both exterior and interior problems. Its diagonalization allowed us to access the eigenvalues $\mu_{mn}$ and eigenfunctions $v_{mn}$ of $\mathcal{M}$, as well as the Steklov eigenfunctions $V_{mn}$. These eigenfunctions inherited the symmetries of the considered domains; in particular, $V_{mn}$ depend on the angle $\phi$ via the factor $e^{im\phi}$, and are symmetric (resp., antisymmetric) with respect to the horizontal plane $z = 0$ when $m+n$ is even (resp., odd). As a consequence, they are also the eigenfunctions of the mixed Steklov-Neumann (resp., Steklov-Dirichlet) problems in the upper half-space. We also described recurrence relations to speed up the numerical construction of truncated matrices. While the interior spectral problem in a bounded domain could alternatively be solved by

other numerical methods [63–73], exterior spectral problems in unbounded domains are in general more difficult to address by conventional tools.

We discussed the impact of domain anisotropy onto the behavior of eigenvalues and eigenfunctions. In particular, we showed how the eigenvalues of $\mathcal{M}$ for the exterior of a prolate spheroid diverge in the limit $a \to 0$ (the exterior of a needle); in turn, in the opposite limit of thin (disk-like) oblate spheroids, the eigenvalues reach finite values. In this limit, we also got complementary insights onto a classical mixed Dirichlet-Neumann problem in the half-space.

Apart from their own fundamental interest in spectral geometry, the Steklov eigenfunctions offer flexible meshless representations for solutions of interior and exterior boundary value problems. In the context of diffusion-controlled reactions, these eigenfunctions allow one to incorporate the effect of partial reactivity, e.g., to compute the steady-state concentration of particles that react on the spheroidal boundary. In particular, we showed how the total diffusive flux (i.e., the overall reaction rate) depends on the anisotropy of the target. Another application from statistical physics concerns the statistics of encounters between a diffusing particle and the spheroidal target. The Steklov eigenfunctions determine the crossing probability $\mathbb{P}_{\boldsymbol{x}_0}\{\ell_\infty > \ell\}$, i.e., the probability that the (rescaled) number of encounters, $\ell_\infty$, exceeds a given threshold $\ell$, before the particle escapes to infinity. The crossing probability generalizes the hitting probability of the target (the latter corresponding to $\ell = 0$). In particular, we showed that the smallest eigenvalue $\mu_{00}$ of the Dirichlet-to-Neumann operator determines the exponential decay of the crossing probability at large thresholds $\ell$.

In summary, this study brings novel insights onto the spectral properties of the Dirichlet-to-Neumann operator for exterior domains, with the special emphasis on anisotropy. Even though the derived eigenvalues and eigenfunctions are not as explicit as in the case of the exterior of a ball, the matrix representation offers an efficient numerical computation and allows for getting asymptotic results. One straightforward extension of this study concerns diffusion across a semi-permeable spheroidal boundary, which is a model of diffusive exchange in anisotropic living cells or tissues [75–77]. The exchange between interior and exterior compartments is usually described by a transmission boundary condition that can be reformulated in terms of two Dirichlet-to-Neumann operators for the interior and exterior problems. Moreover, an extension of the encounter-based approach to semi-permeable membranes allows one to treat much more sophisticated exchange mechanisms [78–81]. Another promising direction is related to the Dirichlet-to-Neumann operator for the (modified) Helmholtz equation. As briefly discussed in Sec. IV B, the related Steklov eigenfunctions determine the encounter propagator and thus most diffusion-reaction characteristics in this system. While one could employ prolate/oblate
<figure><img src="image_14.png" /><figcaption>FIG. 12. Several eigenfunctions $v_{mn}$ of the Dirichlet-to-Neumann operator $M$ for the interior of the prolate spheroid with semi-axes $a = 0.5$ and $b = 1$. The associated eigenvalues are shown on the top. The truncation order is $n_{max} = 10$.</figcaption></figure>

spheroidal wave functions, their analysis is more sophisticated and less explicit. Alternatively, one can consider other domains, in which the Laplace operator admits separation of variables in appropriate curvilinear coordinates (e.g., a torus). A systematic investigation of the Dirichlet-to-Neumann operator in such relatively simple shapes can help to uncover the intricate relations between the geometry and the spectral properties of this operator.

## ACKNOWLEDGMENTS

The author thanks A. Chaigneau for his help on the numerical validation of the proposed method by its comparison to an alternative construction of the Dirichlet-to-Neumann operator via a finite-element method. A partial support by the Alexander von Humboldt Foundation within a Bessel Prize award is acknowledged.

## Appendix A: Interior problems

While the main text is focused on exterior problems, the same approach is valid for the interior spectral problems when the Laplace equation has to solved inside a bounded prolate or oblate spheroid. In fact, the coefficients $F_{n,n'}(z)$ and $\bar{F}_{n,n'}(z)$ remain unchanged, and the difference between exterior and interior settings is only manifested in the coefficients $c_{mn}$:

(i) for prolate spheroids, it is sufficient to replace $Q_n^m(z)$ by $P_n^m(z)$ in Eq. (12), as well as the sign:

$$c_{mn} = \frac{\sinh \alpha_0 P_n^{m'}(\cosh \alpha_0)}{a_E P_n^m(\cosh \alpha_0)} \quad (A1)$$

(ii) for oblate spheroids, it is sufficient to replace $Q_n^m(z)$

<figure><img src="image_15.png" /><figcaption>FIG. 13. Several eigenfunctions $v_{mn}$ of the Dirichlet-to-Neumann operator $M$ for the interior of the oblate spheroid with semi-axes $a = 0.5$ and $b = 1$. The associated eigenvalues are shown on the top. The truncation order is $n_{max} = 10$.</figcaption></figure>

by $P_n^m(z)$ in Eq. (43), as well as the sign:

$$c_{mn} = \frac{i \cosh \alpha_0 P_n^{m'}(i \sinh \alpha_0)}{a_E P_n^m(i \sinh \alpha_0)} \quad (A2)$$

In both cases, as $P_0' = 0$ and thus $c_{00} = 0$, one can easily check that the first eigenvalue is zero, $\mu_{00} = 0$, as it should be.

Figures 12 and 13 illustrate several eigenfunctions $v_{mn}$ for the interior prolate/oblate spheroid with semi-axes $a = 0.5$ and $b = 1$. The associated eigenvalues $\mu_{mn}$ are shown on the top of each panel. Expectedly, a constant eigenfunction $v_{00}$ corresponds to $\mu_{00} = 0$. Other eigenfunctions also resemble the spherical harmonics.

One can also consider more sophisticated domains confined between two confocal prolate (oblate) spheroids. These domains are characterized as $\alpha_1 < \alpha < \alpha_2$, where $\alpha_1 = \tanh^{-1}(a/b)$ and $\alpha_2 = \tanh^{-1}(A/B)$, with $a, b$ and $A, B$ being the semi-axes of the inner and outer boundaries such that $a_E = \sqrt{b^2 - a^2} = \sqrt{B^2 - A^2}$ is half of the focal distance for both boundaries. In this case, a general solution of the Laplace equation involves both $P_n^m(z)$ and $Q_n^m(z)$ but the structure of the matrix representation is very similar. We just mention two settings, in which either Dirichlet or Neumann boundary condition is imposed on the outer boundary, whereas the inner boundary has the Steklov condition (the opposite setting can be easily obtained).

(i) For prolate spheroids, the matrix $M$ is still given by Eq. (14), in which $\alpha_0$ is replaced by $\alpha_1$, and the coefficients $c_{mn}$ depend on the boundary condition on the outer spheroid:
$$c_{mn} = - \frac{\sinh \alpha_1}{a_E} \frac{P_n^m(\cosh \alpha_2) Q_n^{m'}(\cosh \alpha_1) - Q_n^m(\cosh \alpha_2) P_n^{m'}(\cosh \alpha_1)}{P_n^m(\cosh \alpha_2) Q_n^m(\cosh \alpha_1) - Q_n^m(\cosh \alpha_2) P_n^{m'}(\cosh \alpha_1)} \quad (\text{Dirichlet}), \quad (\text{A3})$$

$$c_{mn} = - \frac{\sinh \alpha_1}{a_E} \frac{P_n^{m'}(\cosh \alpha_2) Q_n^{m'}(\cosh \alpha_1) - Q_n^{m'}(\cosh \alpha_2) P_n^{m'}(\cosh \alpha_1)}{P_n^{m'}(\cosh \alpha_2) Q_n^m(\cosh \alpha_1) - Q_n^{m'}(\cosh \alpha_2) P_n^m(\cosh \alpha_1)} \quad (\text{Neumann}). \quad (\text{A4})$$

(ii) For oblate spheroids, the matrix M is still given by Eq. (42), in which $\alpha_0$ is replaced by $\alpha_1$, and the

coefficients $c_{mn}$ depend on the boundary condition on the outer spheroid:

$$c_{mn} = \frac{\cosh \alpha_1}{ia_E} \frac{P_n^m(i \sinh \alpha_2) Q_n^{m'}(i \sinh \alpha_1) - Q_n^m(i \sinh \alpha_2) P_n^{m'}(i \sinh \alpha_1)}{P_n^m(i \sinh \alpha_2) Q_n^m(i \sinh \alpha_1) - Q_n^m(i \sinh \alpha_2) P_n^m(i \sinh \alpha_1)} \quad (\text{Dirichlet}), \quad (\text{A5})$$

$$c_{mn} = \frac{\cosh \alpha_1}{ia_E} \frac{P_n^{m'}(i \sinh \alpha_2) Q_n^m(i \sinh \alpha_1) - Q_n^m(i \sinh \alpha_2) P_n^{m'}(i \sinh \alpha_1)}{P_n^{m'}(i \sinh \alpha_2) Q_n^m(i \sinh \alpha_1) - Q_n^{m'}(i \sinh \alpha_2) P_n^m(i \sinh \alpha_1)} \quad (\text{Neumann}). \quad (\text{A6})$$

## Appendix B: Numerical computation

In this section, we discuss a practical implementation of the proposed method. The construction of the (truncated) matrix M requires computing the integral (15) and the coefficients $c_{mn}$, both involving associated Legendre functions $P_n^m(x)$ and $Q_n^m(x)$. For this reason, we briefly summarize the main steps in the evaluation of these functions and related integrals.

### 1. Associated Legendre functions

The definition and basic properties of associated Legendre functions $P_n^m(x)$ and $Q_n^m(x)$ can be found in many textbooks, e.g., in [82] (chapter 8). As we focus on integer indices $n = 0, 1, 2, ...$ and $m = -n, -n+1, ..., n$, the numerical computation of $P_n^m(x)$ is fairly standard in the common range $-1 \le x \le 1$ (see, e.g., the built-in function `legendre` in matlab). However, the computation of the coefficients $c_{mn}$ requires evaluating both $P_n^m(z)$ and $Q_n^m(z)$ at $z = \cosh \alpha > 1$ or $z = i \sinh \alpha$. While $P_n^m(z)$ with even $m$ are just polynomials and thus can be immediately evaluated at any complex number, the computation is more subtle for other associated Legendre functions that involve square roots and logarithms and thus require cuts in the complex plane. Throughout this section, we follow Ref. [82] and distinguish two conventions by writing the argument as $x$ for $x \in (-1, 1)$, or as $z$ for $z \notin [-1, 1]$. For instance, two conventional definitions read for $m \ge 0$ as:

$$P_n^m(x) = (-1)^m (1 - x^2)^{m/2} \frac{d^m}{dx^m} P_n(x), \quad (\text{B1a})$$

$$P_n^m(z) = (z^2 - 1)^{m/2} \frac{d^m}{dz^m} P_n(z), \quad (\text{B1b})$$

and

$$Q_n^m(x) = (-1)^m (1 - x^2)^{m/2} \frac{d^m}{dx^m} Q_n(x), \quad (\text{B2a})$$

$$Q_n^m(z) = (z^2 - 1)^{m/2} \frac{d^m}{dz^m} Q_n(z), \quad (\text{B2b})$$

where $P_n(z)$ are the Legendre polynomials and $Q_n(z)$ are the Legendre functions of the second kind (see below). In the following, we focus on $P_n^m(z)$ and $Q_n^m(z)$ and briefly describe their numerical computation via the standard recurrence formulas for completeness. As these formulas are identical for both $P_n^m(z)$ and $Q_n^m(z)$, we consider the latter functions and then just mention changes for $P_n^m(z)$.

In the first step, one can compute the Legendre functions $Q_n(z)$ up to the desired order $n$ via the recurrence relation,

$$Q_n(z) = \frac{2n-1}{n} z Q_{n-1}(z) - \frac{n-1}{n} Q_{n-2}(z), \quad (\text{B3})$$

starting from

$$Q_0(z) = \frac{1}{2} \ln \left( \frac{z+1}{z-1} \right), \quad Q_1(z) = zQ_0(z) - 1. \quad (\text{B4})$$

From Eq. (B2b), one can also evaluate

$$Q_n^1(z) = \frac{1}{n\sqrt{z^2-1}} \left( n(2n-1)zQ_n(z) - [(2n-1)(n-1)z^2 + n^2]Q_{n-1}(z) + (n-1)^2zQ_{n-2}(z) \right) \quad (\text{B5})$$

for $n \ge 2$, where we used another recurrence relation

$$\frac{d}{dz} Q_n(z) = \frac{n+1}{z^2-1} (Q_{n+1}(z) - zQ_n(z)). \quad (\text{B6})$$
One also has

$$Q_1^1(z) = \sqrt{z^2 - 1} Q_0(z) - \frac{z}{\sqrt{z^2 - 1}}. \quad (B7)$$

Since we know $Q_n^0(z) = Q_n(z)$ and $Q_n^1(z)$, the other functions $Q_n^m(z)$ can be found via the recurrence relation:

$$Q_n^{m-1}(z) = \frac{2mz(z^2 - 1)^{-\frac{1}{2}} Q_n^m(z) + Q_n^{m+1}(z)}{(n+m)(n-m+1)} \quad (B8)$$

for $m = -n + 2, -n + 3, \dots, -1, 0$, or

$$Q_n^{m+1}(z) = (n+m)(n-m+1)Q_n^{m-1}(z) - \frac{2mzQ_n^m(z)}{(z^2-1)^{\frac{1}{2}}} \quad (B9)$$

for $m = 1, 2, \dots, n$. Note that the derivative of $Q_n^m(z)$ can be found via either of two equivalent relations

$$\frac{dQ_n^m(z)}{dz} = \frac{1}{z^2-1} \left[ (n+m)(n-m+1)\sqrt{z^2-1} Q_n^{m-1}(z) - mzQ_n^m(z) \right], \quad (B10)$$

$$= \frac{1}{z^2-1} \left[ \sqrt{z^2-1} Q_n^{m+1}(z) + mzQ_n^m(z) \right]. \quad (B11)$$

The computation of $P_n^m(z)$ relies on the same recurrence relations, except for their initiation, namely, Eqs. (B4, B7) are replaced by

$$P_0(z) = 1, \quad P_1(z) = z, \quad P_1^1(z) = \sqrt{z^2 - 1}. \quad (B12)$$

## 2. Alternative matrix representation

In the limit $\alpha_0 \to 0$, the coefficients $\bar{F}_{n,n'}^m(\sinh \alpha_0)$ given by Eq. (44) logarithmically diverge that makes their numerical computation more subtle. It is therefore convenient to get an alternative matrix representation that can be constructed even at $\alpha_0 = 0$. For this purpose, we apply the operator $M$ to Eq. (46) that yields

$$\mu_{mn} \sum_{n'=0}^{\infty} [\mathbf{V}_{mn}]_{mn'} \bar{Y}_{mn'} = \mu_{mn} v_{mn} = M v_{mn} \\ = \sum_{n'=0}^{\infty} [\mathbf{V}_{mn}]_{mn'} c_{mn'} \frac{\bar{Y}_{mn'}}{\sqrt{\cosh^2 \alpha_0 - \cos^2 \theta}},$$

where we wrote explicitly the action of $M$ onto a spherical harmonic $\bar{Y}_{mn'}$ via Eq. (41). Multiplying both sides

of this relation by $\sqrt{\cosh^2 \alpha_0 - \cos^2 \theta} \bar{Y}_{mn_0}^*$ cos $\theta$ and integrating over $\theta$ and $\phi$, we get

$$\mu_{mn} \sum_{n'=0}^{\infty} [\mathbf{V}_{mn}]_{mn'} \bar{G}_{mn',mn_0}(\sinh \alpha_0) = [\mathbf{V}_{mn}]_{mn_0} c_{mn_0}, \quad (B13)$$

where we used the orthogonality of spherical harmonics, and defined

$$\begin{aligned} \bar{G}_{mn,m'n'}(z) &= \int_0^{2\pi} d\phi \int_{-\pi/2}^{\pi/2} d\theta \cos\theta \bar{Y}_{mn}(\theta, \phi) \\ & \times \sqrt{z^2 + \sin^2\theta} \bar{Y}_{m'n'}^*(\theta, \phi) \\ &= 2\pi \delta_{m,m'} a_{mn} a_{m'n'} G_{n,n'}^m(z), \end{aligned} \quad (B14)$$

with

$$G_{n,n'}^m(z) = \int_{-1}^{1} dx P_n^m(x) P_{n'}^m(x) \sqrt{z^2 + x^2}. \quad (B15)$$

These coefficients resemble $\bar{F}_{n,n'}^m(z)$ from Eq. (44), except that the square root $\sqrt{z^2 + x^2}$ stands in the numerator and thus eliminates the divergence at $z = 0$. Finally, Eq. (B13) can be written as an eigenvalue problem

$$\sum_{n'=0}^{\infty} [\mathbf{V}_{mn}]_{mn'} \frac{\bar{G}_{mn',mn_0}(\sinh \alpha_0)}{c_{mn_0}} = \frac{1}{\mu_{mn}} [\mathbf{V}_{mn}]_{mn_0}. \quad (B16)$$

This relation implies that $\mathbf{V}_{mn}$ is a left eigenvector of the matrix with the elements $\bar{G}_{mn',mn_0}(\sinh \alpha_0)/c_{mn_0}$, whereas $1/\mu_{mn}$ is the associated eigenvalue. One sees that this is a matrix representation of $M^{-1}$, the inverse of the Dirichlet-to-Neumann operator, which is called the Neumann-to-Dirichlet operator. The numerical advantage of this matrix representation is that the elements $\bar{G}_{n,n'}^m(z)$ do not diverge as $z \to 0$. In particular, one can compute the matrix elements $\bar{G}_{mn',mn_0}(0)/c_{mn_0}$ and then diagonalize the truncated matrix.

## 3. Recurrence formulas for the integral

Now we describe the recurrence formulas for computing the elements of the matrix $M$. These formulas are actually valid in a more general case:

$$H_{n,n'}^m = \int_{-1}^{1} dx f(x^2) P_n^m(x) P_{n'}^m(x), \quad (B17)$$

where $f(x)$ is a given integrable function. Due to the symmetry relation (25), these integrals vanish when $n+n'$ is odd. They also are zero if $|m| > n$ or $|m| > n'$. The remaining nonzero values can be found recursively by using the recurrence relation for associated Legendre polynomials:

$$(n-m+1)P_{n+1}^m(x) = (2n+1)xP_n^m(x) - (n+m)P_{n-1}^m(x). \quad (B18)$$

Applying this relation to both $P_n^m(x)$ and $P_{n'}^m(x)$ in Eq. (B17), we get for even $n+n'$:
$$
\begin{align}
H_{n,n'}^m &= \int_{-1}^{1} dx f(x^2) P_n^m(x) \frac{(2n' - 1)x P_{n'-1}^m(x) - (n' - 1 + m) P_{n'-2}^m(x)}{n' - m} \nonumber \\
&= \int_{-1}^{1} dx f(x^2) \left[ \frac{(n - m + 1)P_{n+1}^m(x) + (n + m)P_{n-1}^m(x)}{2n + 1} (2n' - 1)P_{n'-1}^m(x) - \frac{n' - 1 + m}{n' - m} P_n^m(x) P_{n'-2}^m \right] \nonumber \\
&= \frac{2n' - 1}{(2n + 1)(n' - m)} \left( (n + 1 - m)H_{n+1,n'-1}^m + (n + m)H_{n-1,n'-1}^m \right) - \frac{n' - 1 + m}{n' - m} H_{n,n'-2}^m. \tag{B19}
\end{align}
$$

As a consequence, if $H_{n,0}^m$ and $H_{n,1}^m$ are known, one can
evaluate the remaining elements.

For any $m > 0$, one has to evaluate

$$
\begin{align*}
H_{n,m}^m &= 2 \int_0^1 dx f(x^2) P_n^m(x) \underbrace{(-1)^m (2m-1)!! (1-x^2)^{m/2}}_{=P_m^m(x)} \\
&= -(2m-1)2 \int_0^1 dx f(x^2) P_{m-1}^{m-1}(x) \sqrt{1-x^2} P_n^m(x) \\
&= \frac{2m-1}{2n+1} \left( (n+m-1)(n+m)H_{n-1,m-1}^{m-1} \right. \\
&\qquad \left. -(n-m+1)(n-m+2)H_{n+1,m-1}^{m-1} \right)
\end{align*}
$$

for even $m + n$, where we used another recurrence relation to express $\sqrt{1-x^2}P_n^m(x)$ in terms of $P_{n\pm 1}^{m-1}(x)$. Similarly, one needs

$$
H_{n,m+1}^m = \frac{2m+1}{2n+1} \left( (n+m-1)(n+m)H_{n-1,m}^{m-1} - (n-m+1)(n-m+2)H_{n+1,m}^{m-1} \right)
$$

for odd $m+n$. In other words, once the elements $H_{n,n'}^{m-1}$ are constructed, one can find the elements $H_{n,n'}^m$. As a consequence, if one knows the elements $H_{n,0}^0$ and $H_{n,1}^0$, one can first construct all the elements $H_{n,n'}^0$ by using recurrence relations (B19), and then progressively get the elements $H_{n,n'}^m$. We stress that this procedure does not depend on the function $f$, which determines only the initialization step, i.e., the elements $H_{n,0}^0$ and $H_{n,1}^0$.

For the integrals $\bar{F}_{n,n'}^m(z)$ from Eq. (44), one has
$f(y) = (z^2 + y)^{-1/2}$, for which we have

$$
\bar{F}_{0,0}^{0}(z) = \ln \left( \frac{\sqrt{z^2 + 1} + 1}{\sqrt{z^2 + 1} - 1} \right), \qquad (B20a)
$$

$$
\bar{F}_{1,1}^{0}(z) = \sqrt{z^2 + 1} - \frac{z^2}{2} \ln \left( \frac{\sqrt{z^2 + 1} + 1}{\sqrt{z^2 + 1} - 1} \right). \quad (B20b)
$$

Using the recurrence relations

$$
P_n(x) = \frac{2n-1}{n} x P_{n-1}(x) - \frac{n-1}{n} P_{n-2}(x) \quad (B21)
$$

for Legendre polynomials $P_n(x)$ and integrating by parts,
one gets two sets of relations

$$
\bar{F}_{n+1,0}^{0}(z) = \frac{2n+1}{n+1} \bar{F}_{n,1}^{0}(z) - \frac{n}{n+1} \bar{F}_{n-1,0}^{0}(z), \quad (B22a)
$$

$$
\bar{F}_{n+1,1}^{0}(z) = -\frac{n-1}{n+2} \bar{F}_{n-1,1}^{0}(z) - \frac{2n+1}{n+2} z^2 \bar{F}_{n,0}^{0}(z). \quad (B22b)
$$

These relations determine all $\bar{F}_{n,0}^0(z)$ and $\bar{F}_{n,1}^0(z)$, from
which one can construct recursively all $\bar{F}_{n,n'}^m(z)$, as de-
scribed above. According to Eq. (45), this construction
can also be applied to get the elements $F_{n,n'}^m(z)$.

For the integrals $\tilde{G}_{n,n'}^m(z)$ from Eq. (B15), one has
$f(y) = (z^2 + y)^{1/2}$, for which we have

$$
\bar{G}_{0,0}^{0}(z) = \sqrt{z^2 + 1} + \frac{z^2}{2} \ln \left( \frac{\sqrt{z^2 + 1} + 1}{\sqrt{z^2 + 1} - 1} \right), \quad (B23a)
$$

$$
\bar{G}_{1,1}^{0}(z) = \frac{(z^2 + 2)\sqrt{z^2 + 1}}{4} - \frac{z^4}{8} \ln \left( \frac{\sqrt{z^2 + 1} + 1}{\sqrt{z^2 + 1} - 1} \right). \quad (B23b)
$$

Using the recurrence relations and integrating by parts,
one retrieves Eq. (B22a), which is completed by

$$
\bar{G}_{n+1,1}^{0}(z) = -\frac{n-3}{n+4} \bar{G}_{n-1,1}^{0}(z) - \frac{2n+1}{n+4} z^2 \bar{G}_{n,0}^{0}(z). \quad (B24)
$$

These relations determine all $\bar{G}_{n,0}^0(z)$ and $\bar{G}_{n,1}^0(z)$, from which one can construct recursively all $\bar{G}_{n,n'}^m(z)$, as described above.

In the limit of a disk ($z = \sinh \alpha_0 \to 0$), one gets
$\bar{G}_{0,0}^0(0) = 1$ and $\bar{G}_{1,1}^0(0) = 1/2$, from which all other
elements $\bar{G}_{n,n'}^m(0)$ can be recursively constructed.

4. Orthogonality of eigenfunctions

Let us check that the eigenfunctions $v_k$ obtained from
the eigenvectors $\mathbf{V}_k$ of the matrix $\mathbf{M}$ are orthogonal to
each other. For oblate spheroids, we have

$$
\begin{align}
(v_{k_1}, v_{k_2})_{L^2(\partial\Omega)} &= \int_{-\pi/2}^{\pi/2} d\theta \int_0^{2\pi} d\phi h_\theta h_\phi v_{k_1}^* v_{k_2} \nonumber \\
&= a_E^2 \cosh \alpha_0 \sum_{m_1, n_1, m_2, n_2} [\mathbf{V}_{k_1}^*]_{m_1 n_1} [\mathbf{V}_{k_2}]_{m_2 n_2} \nonumber \\
&\quad \times \bar{G}_{m_1 n_1, m_2 n_2} (\sinh \alpha_0) \nonumber \\
&= a_E^2 \cosh \alpha_0 \frac{1}{\mu_{k_1}} \sum_{m_2, n_2} [\mathbf{V}_{k_1}^*]_{m_2 n_2} c_{m_2 n_2} [\mathbf{V}_{k_2}]_{m_2 n_2}, \tag{B25}
\end{align}
$$

where we used Eqs. (B13, B16), and $\bar{G}_{m_1 n_1, m_2 n_2} (\sinh \alpha_0)$ are given by Eq. (B14). It remains to show that the above sum vanishes when $\mu_{k_1} \neq \mu_{k_2}$.

For this purpose, we rewrite Eq. (42) in a matrix form as $\mathbf{M} = \mathbf{F}\mathbf{c}$, where $\mathbf{F}$ is the symmetric matrix with elements

$$
\bar{F}_{mn, m'n'} (\sinh \alpha_0) = 2\pi \delta_{m,m'} a_{mn} a_{m'n'} \bar{F}_{n,n'}^m (\sinh \alpha_0), \quad (B26)
$$

and $\mathbf{c}$ is the diagonal matrix formed by $c_{mn}$. On one hand, one can transpose the relation $\mathbf{M}\mathbf{V}_{k_1} = \mu_{k_1}\mathbf{V}_{k_1}$ and multiply it by $\mathbf{c}\mathbf{V}_{k_2}$ on the right to get

$$
\mathbf{V}_{k_1}^{\dagger} \mathbf{c} \mathbf{F} \mathbf{c} \mathbf{V}_{k_2} = \mu_{k_1} \mathbf{V}_{k_1}^{\dagger} \mathbf{c} \mathbf{V}_{k_2}.
$$

On the other hand, multiplication of the relation
$\mathbf{M}\mathbf{V}_{k_2} = \mu_{k_2}\mathbf{V}_{k_2}$ by $\mathbf{V}_{k_1}^{\dagger}\mathbf{c}$ on the left yields

$$
\mathbf{V}_{k_1}^{\dagger} \mathbf{c} \mathbf{F} \mathbf{c} \mathbf{V}_{k_2} = \mu_{k_2} \mathbf{V}_{k_1}^{\dagger} \mathbf{c} \mathbf{V}_{k_2}.
$$

Subtracting these equations, one gets

$$
(\mu_{k_1} - \mu_{k_2}) \mathbf{V}_{k_1}^{\dagger} \mathbf{c} \mathbf{V}_{k_2} = 0,
$$

that implies the orthogonality for any $\mu_{k_1} \neq \mu_{k_2}$:

$$
\mathbf{V}_{k_1}^{\dagger} \mathbf{c} \mathbf{V}_{k_2} = \sum_{m,n} [\mathbf{V}_{k_1}^{*}]_{mn} c_{mn} [\mathbf{V}_{k_2}]_{mn} = 0, \quad (\text{B27})
$$

and thus the orthogonality of the eigenfunctions $v_{k_1}$ and
$v_{k_2}$ due to Eq. (B25).

Since the matrix **M** is not symmetric, its eigenvectors
**V**<sub>k</sub> are not orthogonal to each other; in fact, their orthog-
onality relation (B27) includes the weighting coefficient
**c**<sub>mn</sub>. For this reason, it is more convenient to consider
rescaled eigenvectors **V**̃<sub>k</sub> = **c**<sub>1</sub><sup>2</sup> **V**<sub>k</sub>, which are the eigen-
 vectors of the Hermitian matrix

$$
\tilde{\mathbf{M}} = \mathbf{c}^{\frac{1}{2}} \mathbf{M} \mathbf{c}^{-\frac{1}{2}} = \mathbf{c}^{\frac{1}{2}} \mathbf{F} \mathbf{c}^{\frac{1}{2}}. \quad (\text{B28})
$$

In fact, one has $\tilde{\mathbf{M}}\tilde{\mathbf{V}}_k = \mu_k \tilde{\mathbf{V}}_{k_z}$ with the same eigenvalue $\mu_k$. The rescaled eigenvectors $\tilde{\mathbf{V}}_k$ are orthogonal to each.

Note that Eq. (B25) allows one to compute the
$L^2(\partial\Omega)$-norm of each eigenfunction directly from the cor-

responding eigenvector:

$$
\begin{align*}
\|v_k\|_{L^2(\partial\Omega)}^2 &= \int_{-\pi/2}^{\pi/2} d\theta \int_0^{2\pi} d\phi h_\theta h_\phi |v_k|^2 \\
&= a_E^2 \cosh \alpha_0 \frac{1}{\mu_k} \sum_{m,n} [\mathbf{V}_k^*]_{mn} c_{mn} [\mathbf{V}_k]_{mn} \\
&= a_E^2 \cosh \alpha_0 \frac{\mathbf{V}_k^\dagger \mathbf{c} \mathbf{V}_k}{\mu_k} = a_E^2 \cosh \alpha_0 \frac{\tilde{\mathbf{V}}_k^\dagger \tilde{\mathbf{V}}_k}{\mu_k}. \tag{B29}
\end{align*}
$$

This straightforward computation helps to avoid numeri-
cal integration over θ. A similar computation for prolate
spheroids yields:

$$
\|v_k\|_{L^2(\partial\Omega)}^2 = a_E^2 \sinh \alpha_0 \frac{\mathbf{V}_k^\dagger \mathbf{c} \mathbf{V}_k}{\mu_k} = a_E^2 \sinh \alpha_0 \frac{\tilde{\mathbf{V}}_k^\dagger \tilde{\mathbf{V}}_k}{\mu_k}. \quad (B30)
$$

In the same vein, one can compute the projection of an
eigenfunction $v_k$ onto a constant. For oblate spheroids,
one gets

$$
\begin{align*}
(1, v_k)_{L^2(\partial\Omega)} &= \int_{-\pi/2}^{\pi/2} d\theta \int_0^{2\pi} d\phi h_\theta h_\phi v_k \\
&= a_E^2 \cosh \alpha_0 \sum_{m,n} [\mathbf{V}_k]_{mn} \bar{G}_{mn,00} (\sinh \alpha_0) \sqrt{4\pi} \\
&= \sqrt{4\pi} a_E^2 \cosh \alpha_0 \frac{c_{00}}{\mu_k} [\mathbf{V}_k]_{00}, \tag{B31}
\end{align*}
$$

where we used that $\tilde{Y}_{00} = a_{00} = 1/\sqrt{4\pi}$. For instance, for
the exterior problem, substitution of $c_{00}$ from Eq. (43)
yields

$$
(1, v_k)_{L^2(\partial\Omega)} = \frac{\sqrt{4\pi a_E}}{iQ_0(i \sinh \alpha_0)} \frac{[\mathbf{V}_k]_{00}}{\mu_k}. \quad (B32)
$$

For prolate spheroids, a similar computation gives

$$
(1, v_k)_{L^2(\partial\Omega)} = \sqrt{4\pi} a_E^2 \sinh \alpha_0 \frac{c_{00}}{\mu_k} [\mathbf{V}_k]_{00}. \quad (B33)
$$

For instance, for the exterior problem, substitution of $c_{00}$ from Eq. (12) yields

$$
(1, v_k)_{L^2(\partial\Omega)} = \frac{\sqrt{4\pi} a_E}{Q_0(\cosh \alpha_0)} \frac{[\mathbf{V}_k]_{00}}{\mu_k}. \quad (B34)
$$

**5. Axisymmetric setting**

In many applications, it is sufficient to look at ax-
isymmetric eigenfunctions, which do not depend on the
angle φ and thus correspond to m = 0. These eigen-
functions can be constructed by diagonalizing the matrix
M₀, for which the computation can be further simplified.
We briefly discuss this case below for prolate and oblate
spheroids.
Prolate spheroids

In this case, Eq. (10) is reduced to

$$f(\theta) = \sum_{n=0}^{\infty} f_n \psi_n(\theta), \qquad (B35)$$

where

$$\psi_n(\theta) = \sqrt{n+1/2} P_n(\cos \theta) \qquad (B36)$$

are normalized Legendre polynomials. The action of $M$ reads then

$$[Mf](\theta) = \sum_{n=0}^{\infty} [M_0 f]_n \psi_n(\theta), \qquad (B37)$$

where the (infinite-dimensional) matrix $M_0$ represents the operator $M$ on the orthonormal basis of Legendre polynomials, with

$$[M_0]_{n,n'} = M_{0,n,0n'} = \sqrt{(n+1/2)(n'+1/2)} c_{0n'} F_{n,n'}^0(\cosh \alpha_0), \qquad (B38)$$

where

$$c_{0n} = - \frac{\sinh \alpha_0 Q'_n (\cosh \alpha_0)}{a_E Q_n (\cosh \alpha_0)}, \qquad (B39)$$

and

$$F_{n,n'}^0(z) = \int_{-1}^{1} dx \frac{P_n(x) P_{n'}(x)}{\sqrt{z^2 - x^2}}. \qquad (B40)$$

By symmetry, the elements $F_{n,n'}^0(z)$ are zero when $n+n'$ is odd. Once the eigenvectors $V_{0n}$ of the matrix $M_0$ are

found, the eigenfunction $v_{0n}$ and the Steklov eigenfunction $V_{0n}$ are given by

$$v_{0n}(\theta) = \sum_{n'=0}^{\infty} [V_{0n}]_{0n'} \psi_{n'}(\theta) \qquad (B41)$$

and

$$V_{0n}(\alpha, \theta) = \sum_{n'=0}^{\infty} \frac{Q_{n'}(\cosh \alpha)}{Q_{n'}(\cosh \alpha_0)} [V_{0n}]_{0n'} \psi_{n'}(\theta). \qquad (B42)$$

Oblate spheroids

In the same vein, one employs the normalized Legendre polynomials

$$\psi_n(\theta) = \sqrt{n+1/2} P_n(\sin \theta) \qquad (B43)$$

to get the matrix representation of $M$:

$$[M_0]_{0n,0n'} = \sqrt{(n+1/2)(n'+1/2)} c_{0n'} i F_{n,n'}^0(i \sinh \alpha_0), \qquad (B44)$$

where

$$c_{0n} = \frac{\cosh \alpha_0 Q'_n (i \sinh \alpha_0)}{i a_E Q_n (i \sinh \alpha_0)}. \qquad (B45)$$

Disk

In the limit $a=0$, the coefficients $c_{0n}$ from Eq. (B45) approach constants; in fact, using

$$Q_{2n}(0) = (-1)^n \frac{(2n)!}{4^n (n!)^2} \frac{\pi}{2i}, \qquad (B46)$$

$$Q_{2n-1}(0) = (-1)^n \frac{2^{2n-1} (n!)^2}{n(2n)!}, \qquad (B47)$$

and the recurrence relation (B6) to evaluate $Q'_n(0) = nQ_{n-1}(0)$, we deduce Eq. (51). Note that these coefficients can be found iteratively as

$$c_{0n} = \frac{c_{0(n-2)}}{(1-1/n)^2}. \qquad (B48)$$

One has $c_{0n} \simeq n$ as $n$ increases.

[1] M. Cheney, D. Isaacson, and J. C. Newell, Electrical Impedance Tomography, SIAM Rev. **41**, 85-101 (1999).

[2] A. P. Calderón, On an inverse boundary value problem, Seminar on Numerical Analysis and its Applications to Continuum Physics, Soc. Brasileira de Matemática, Río de Janeiro, (1980), 65-73; Reprinted in Comput. Appl. Math. **25**, 2-3, 133-138 (2006).

[3] L. Borcea, Electrical impedance tomography, Inv. Prob. **18**, R99 (2002).

[4] J. Sylvester and G. Uhlmann, A Global Uniqueness Theorem for an Inverse Boundary Value Problem, Ann. Math. **125**, 153-169 (1987).

[5] E. Curtis and J. Morrow, The Dirichlet to Neumann map for a resistor network, SIAM J. Appl. Math. **51**, 1011-1029 (1991).

[6] M. S. Agranovich, B. Z. Katsenelenbaum, A. N. Sivov, and N. N. Voitovich, *Generalized Method of Eigenoscillations in Diffraction Theory* (Wiley-VCH, Berlin, 1999).
[7] B. F. Smith, Domain Decomposition Methods for Partial Differential Equations, in "Parallel Numerical Algorithms", Eds. D. E. Keyes, A. Sameh, V. Venkatakrishnan (Springer, 1996), pp. 225-243.

[8] D. Givoli, I. Patlashenko, and J. B. Keller, Discrete Dirichlet-to-Neumann maps for unbounded domains, Comput. Methods Appl. Mech. Engrg. **164**, 173-185 (1998).

[9] A. Delitsyn and D. S. Grebenkov, Mode matching methods in spectral and scattering problems, Quart. J. Mech. Appl. Math. **71**, 537-580 (2018).

[10] A. Delitsyn and D. S. Grebenkov, Resonance scattering in a waveguide with identical thick perforated barriers, Appl. Math. Comput. **412**, 126592 (2022).

[11] G. Auchmuty, Steklov eigenproblems and the representation of solutions of elliptic boundary value problems, Numer. Funct. Anal. Optim. **25**, 321-348 (2004).

[12] G. Auchmuty and Q. Han, Spectral representations of solutions of linear elliptic equations on exterior regions, J. Math. Anal. Appl. **398**, 1-10 (2013).

[13] G. Auchmuty and Q. Han, Representations of Solutions of Laplacian Boundary Value Problems on Exterior Regions, Appl. Math. Optim. **69**, 21-45 (2014).

[14] G. Auchmuty and M. Cho, Boundary integrals and approximations of harmonic functions, Numer. Funct. Anal. Optim. **36**, 687-703 (2015).

[15] G. Auchmuty, Steklov Representations of Green's Functions for Laplacian Boundary Value Problems, Appl. Math. Optim. **77** (2018).

[16] D. S. Grebenkov, M. Filoche, and B. Sapoval, Mathematical Basis for a General Theory of Laplacian Transport towards Irregular Interfaces, Phys. Rev. E **73**, 021103 (2006).

[17] D. S. Grebenkov, Spectral theory of imperfect diffusion-controlled reactions on heterogeneous catalytic surfaces, J. Chem. Phys. **151**, 104108 (2019).

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

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

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

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

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

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

[24] M. Levitin, D. Mangoubi, and I. Polterovich, Topics in Spectral Geometry (Graduate Studies in Mathematics, vol. 237; American Mathematical Society, 2023).

[25] A. Girouard and I. Polterovich, Spectral geometry of the Steklov problem, J. Spectr. Th. **7**, 321-359 (2017).

[26] B. Colbois, A. Girouard, C. Gordon, and D. Sher, Some recent developments on the Steklov eigenvalue problem, Rev. Mat. Complut. (2023).

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

[28] W. Arendt and A. F. M. ter Elst, The Dirichlet-to-Neumann Operator on Exterior Domains, Potential Anal. **43**, 313-340 (2015).

[29] T. J. Christiansen and K. Datchev, Low energy scattering asymptotics for planar obstacles, Pure Appl. Anal. **5**, 767-794 (2023).

[30] C. Xiong, Sharp bounds for the first two eigenvalues of an exterior Steklov eigenvalue problem, ArXiv:2304.11297v1 (2023).

[31] P. M. Morse and H. Feshbach, Methods of Theoretical Physics (New York: McGraw-Hill, 1953).

[32] W. R. Smythe, *Static and Dynamic Electricity*, 3rd ed. (New York: McGraw-Hill, 1968).

[33] H. Barucq, R. Djellouli, and A. Saint-Guirons, Three-dimensional approximate local DtN boundary conditions for prolate spheroid boundaries, J. Comput. Appl. Math. **234**, 1810-1816 (2010).

[34] Q. T. Le Gia, E. P. Stephan, and T. Tran, Solution to the Neumann problem exterior to a prolate spheroid by radial basis functions, Adv. Comput. Math. **34**, 83-103 (2011).

[35] A. Costea, Q. T. Le Gia, E. P. Stephan, and T. Tran, Meshless BEM and overlapping Schwarz preconditioners for exterior problems on spheroids, Stud. Geophys. Geod. **55**, 465-477 (2011).

[36] X. Luo, Q. Du, and L. Liu, A D-N alternating algorithm for exterior 3D Poisson problem with prolate spheroid boundary, Appl. Math. Comput. **269**, 252-264 (2015).

[37] D. Gomez and A. V. Cheviakov, Asymptotic analysis of narrow escape problems in nonspherical three-dimensional domains, Phys. Rev. E **91**, 012137 (2015).

[38] C. Xue ans S. Deng, Green's function and image system for the Laplace operator in the prolate spheroidal geometry, AIP Advances **7**, 015024 (2017).

[39] S. D. Traytak and D. S. Grebenkov, Diffusion-influenced reaction rates for active "sphere-prolate spheroid" pairs and Janus dimers, J. Chem. Phys. **148**, 024107 (2018).

[40] F. Piazza and D. S. Grebenkov, Diffusion-controlled reaction rate on non-spherical partially absorbing axisymmetric surfaces, Phys. Chem. Chem. Phys. **21**, 25896-25906 (2019).

[41] A. Chaigneau and D. S. Grebenkov, First-passage times to anisotropic partially reactive targets, Phys. Rev. E **105**, 054146 (2022).

[42] A. Chaigneau and D. S. Grebenkov, Effects of target anisotropy on harmonic measure and mean first-passage time, J. Phys. A: Math. Theor. **56**, 235202 (2023).

[43] G. A. Korn and T. M. Korn, *Mathematical Handbook for Scientists and Engineers* (New York: McGraw-Hill, 1961).

[44] I. N. Sneddon, *Mixed Boundary Value Problems in Potential Theory* (Wiley, NY, 1966).

[45] J. B. Garnett and D. E. Marshall, *Harmonic Measure* (Cambridge University Press, 2005).

[46] N. G. Van Kampen, *Stochastic Processes in Physics and Chemistry* (Elsevier, Amsterdam, 1992).

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

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

[49] R. Metzler, G. Oshanin, and S. Redner (Eds), *First-Passage Phenomena and Their Applications* (Singapore, World Scientific, 2014).
[50] K. Lindenberg, R. Metzler, and G. Oshanin (Eds), *Chemical Kinetics: Beyond the Textbook* (World Scientific, New Jersey, 2019).

[51] A. M. North, Diffusion-controlled reactions, Q. Rev. Chem. Soc. **20**, 421–440 (1966).

[52] G. Wilemski and M. Fixman, General theory of diffusion-controlled reactions, J. Chem. Phys. **58**, 4009–4019 (1973).

[53] D. F. Calef and J. M. Deutch, Diffusion-Controlled Reactions, Ann. Rev. Phys. Chem. **34**, 493–524 (1983).

[54] O. G. Berg and P. H. von Hippel, Diffusion-Controlled Macromolecular Interactions, Ann. Rev. Biophys. Bio-phys. Chem. **14**, 131–160 (1985).

[55] S. A. Rice, *Diffusion-limited reactions* (Elsevier, Amsterdam, 1985).

[56] P. C. Bressloff and J. Newby, Stochastic models of intracellular transport, Rev. Mod. Phys. **85**, 135-196 (2013).

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

[58] D. S. Grebenkov, Diffusion-Controlled Reactions: An Overview, Molecules **28**, 7570 (2023).

[59] M. Smoluchowski, Versuch einer Mathematischen Theorie der Koagulations Kinetic Kolloider Lösungen, Z. Phys. Chem. **92U**, 129-168 (1918).

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

[61] D. S. Grebenkov and D. Krapf, Steady-state reaction rate of diffusion-controlled reactions in sheets, J. Chem. Phys **149**, 064117 (2018).

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

[63] A. B. Andreev and T. D. Todorov, Isoparametric finite-element approximation of a Steklov eigenvalue problem, IMA J. Num. Anal. **24**, 309-322 (2004).

[64] Y. Yang, Q. Li, and S. Li, Nonconforming finite element approximations of the Steklov eigenvalue problem, Appl. Num. Math. **59**, 2388-2401 (2009).

[65] H. Bi, Y. Yang, A two-grid method of the non-conforming Crouzeix-Raviart element for the Steklov eigenvalue problem, Appl. Math. Comput. **217**, 9669-9678 (2011).

[66] Q. Li and Y. Yang, A two-grid discretization scheme for the Steklov eigenvalue problem, J. Appl. Math. Comput. **36**, 129-139 (2011).

[67] Q. Li, Q. Lin, and H. Xie, Nonconforming finite element approximations of the Steklov eigenvalue problem and its lower bound approximations, Appl. Math. **58**, 129-151 (2013).

[68] H. Xie, A type of multilevel method for the Steklov eigenvalue problem, IMA J. Num. Anal. **34**, 592-608 (2014).

[69] H. Bi, H. Li, and Y. Yang, An adaptive algorithm based on the shifted inverse iteration for the Steklov eigenvalue problem, Appl. Num. Math. **105** 64-81 (2016).

[70] B. Bogosel, The method of fundamental solutions applied to boundary eigenvalue problems, J. Comput. Appl. Math. **306**, 265-285 (2016).

[71] E. Akhmetgaliyev, C.-Y. Kao, and B. Osting, Computational methods for extremal Steklov problems, SIAM J. Contr. Optim. **55**, 1226-1240 (2017).

[72] W. Alhejaili and C.-Y. Kao, Numerical studies of the Steklov eigenvalue problem via conformal mappings, Appl. Math. Comput. **347**, 785-802 (2019).

[73] J.-T. Chen, J.-W. Lee, and K.-T. Lien, Analytical and numerical studies for solving Steklov eigenproblems by using the boundary integral equation method/boundary element method, Engnr. Anal. Bound. Elem. **114**, 136-147 (2020).

[74] A. Chaigneau and D. S. Grebenkov, A numerical study of the Dirichlet-to-Neumann operator in planar domains (submitted).

[75] P. H. Richter and M. Eigen, Diffusion controlled reaction rates in spheroidal geometry: Application to repressor-operator association and membrane bound enzymes, Bio-phys. Chem. **2**, 255-263 (1974).

[76] D. R. Grimes and F. J. Currell, Oxygen diffusion in ellipsoidal tumour spheroids, J. R. Soc. Interface **15**, 20180256 (2018).

[77] M. Gadzinowski, D. Mickiewicz, and T. Basinska, Spherical versus prolate spheroidal particles in biosciences: Does the shape make a difference? Polym. Adv. Technol. 1-10 (2021).

[78] P. C. Bressloff, Diffusion-mediated absorption by partially-reactive targets: Brownian functionals and generalized propagators, J. Phys. A: Math. Theor. **55**, 205001 (2022).

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

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

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

[82] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables* (United States Department of Commerce, National Bureau of Standards, 1964).