# The Steklov problem for exterior domains: Asymptotic behavior and applications

Denis S Grebenkov, A. Chaigneau

► To cite this version:

Denis S Grebenkov, A. Chaigneau. The Steklov problem for exterior domains: Asymptotic behavior and applications. Journal of Mathematical Physics, 2025, 66 (6), pp.061502. 10.1063/5.0228529. hal-05344136

HAL Id: hal-05344136

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

Submitted on 3 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# The Steklov problem for exterior domains: asymptotic behavior and applications

D. S. Grebenkov¹,² and A. Chaigneau¹

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

²)CNRS – Université de Montréal CRM – CNRS, 6128 succ Centre-Ville, Montréal QC H3C 3J7, Canada

(*Electronic mail: denis.grebenkov@polytechnique.edu)

(Dated: 25 March 2025)

We investigate the spectral properties of the Steklov problem for the modified Helmholtz equation $(p - \Delta)u = 0$ in the exterior of a compact set, for which the positive parameter $p$ ensures exponential decay of the Steklov eigenfunctions at infinity. We obtain the small-$p$ asymptotic behavior of the eigenvalues and eigenfunctions and discuss their features for different space dimensions. These results find immediate applications to the theory of stochastic processes and unveil the long-time asymptotic behavior of probability densities of various first-passage times in exterior domains. Theoretical results are validated by solving the exterior Steklov problem by a finite-element method with a transparent boundary condition.

# I. INTRODUCTION

The Laplace operator plays the central role in mathematics, physics, and many other disciplines¹⁻¹². So, its eigenmodes, satisfying $-\Delta u = \lambda u$ in a considered bounded domain $\Omega \subset \mathbb{R}^d$ with an appropriate boundary condition, can represent vibration modes of a drum, quantum states of a free particle in a confinement (e.g., a waveguide), or proper basis functions for diffusive processes, whereas the Laplacian eigenvalues are respectively related to frequencies, energies, and decay times of these eigenmodes. In a similar spectral problem formulated by Steklov¹³,¹⁴, the spectral parameter stands in the Robin-type boundary condition, i.e., one searches for eigenpairs {$\mu, u$} satisfying

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

where $\partial_n$ is the normal derivative on the boundary $\partial\Omega$ oriented outwards the domain $\Omega$. Despite intensive mathematical studies over the past century (see the recent book¹⁵ and reviews¹⁶,¹⁷), the Steklov problem is much less known in physics, while its physical applications, such as electrical impedance tomography¹⁸⁻²² or sloshing problem²³⁻²⁵, are less numerous. Recently, the eigenmodes of the Steklov problem have found a new application within the encounter-based approach to diffusion-controlled reactions²⁶⁻²⁹. In particular, the stochastic formulation of such diffusion-mediated phenomena has naturally led to spectral expansions related to the Steklov problem.

One of the crucial advantages of the Steklov problem manifests for unbounded Euclidean domains with a bounded boundary. While former spectral expansions based on Laplacian eigenfunctions become unavailable due to the continuum spectrum of the Laplace operator, the spectrum of the Steklov problem is still discrete, allowing one to employ its spectral expansions. For instance, Auchmuty suggested using the Steklov eigenfunctions for representing harmonic functions in exterior domains and for solving associated boundary value problems³⁰⁻³⁴. In diffusion applications, one often deals with particles diffusing in the unbounded space towards bounded targets³⁵⁻⁴⁰ so that the Steklov eigenmodes can be particularly valuable (see, e.g., the statistics of encounters between two diffusing disks⁴¹).

Despite this growing interest, most former works on the Steklov problem were restricted to bounded domains. From the mathematical point of view, the analysis of the Steklov problem is more difficult for exterior problems. For instance, a constant function, which is the principal Steklov eigenmode for any bounded domain, is not square integrable in the case of unbounded domains. More generally, as formal solutions of the exterior Steklov problem may not be square integrable, one needs either to employ larger functional spaces for a proper formulation of the Steklov problem³³,³⁴,⁴², or to consider the interior Steklov problem in a sequence of bounded domains approaching the unbounded one⁴³. From the numerical point of view, one cannot directly apply common tools such as a finite-element method (FEM) to an unbounded domain that requires a truncation of the original domain and dealing with the resulting artificial boundary (see below). From the applicative point of view, the limited knowledge on the exterior Steklov problem delays its versatile applications in physics and other disciplines.

In this paper, we follow another approach to formulate and investigate the exterior Steklov problem. We start from a generalized Steklov problem, in which the Laplace equation is replaced by the modified Helmholtz equation, $(p - \Delta)u = 0$, with a fixed parameter $p > 0$. An exponential decay of solutions of this equation at infinity allows one to handle the generalized Steklov problem in exterior domains in the same conventional way, as for bounded domains³¹ (see also Appendix A for exponentially decaying upper bounds). This spectral problem presents its own mathematical interest and has direct applications to stochastic processes, first-passage phenomena, and diffusion-controlled reactions, as illustrated in Sec. VI. We focus on the asymptotic
behavior of the eigenvalues and eigenfunctions of the generalized exterior Steklov problem in the limit $p \to 0$. In particular, we show how this behavior strongly depends on the dimension $d$ of the embedding Euclidean space $\mathbb{R}^d$. The limits of the eigenvalues and eigenfunctions at $p = 0$ are then *postulated* to be the eigenmodes of the exterior Steklov problem (1). The existence of these limits and the equivalence of this approach with other formulations$^{33,43}$ of the exterior Steklov problem are proved in a separate work$^{44}$.

The paper is organized as follows. In Sec. II, we formulate the generalized exterior Steklov problem and recall its basic mathematical properties. In Sec. III, we inspect the exterior of a ball as a guiding example, for which this spectral problem admits an explicit solution. Its spectral properties bring some intuition and practical tools for studying the general case in Sec. IV. Section V presents some numerical illustrations for two-dimensional and three-dimensional domains. In Sec. VI, we discuss an application of the derived asymptotic results to diffusion-controlled reactions and the related first-passage time statistics. In particular, we derive the long-time behavior of the survival probability for exterior domains in three dimensions. Section VII concludes the paper.

## II. GENERALIZED EXTERIOR STEKLOV PROBLEM

Let $\Omega \in \mathbb{R}^d$ ($d \ge 2$) be the exterior of a compact set, and $\partial\Omega$ denote its smooth connected boundary (while the interior Steklov problem is often formulated for domains with Lipschitz boundary$^{15}$, we deliberately avoid the related subtle points by focusing on smooth boundaries; at the same time, numerical examples in Sec. V include polygonal domains; see also more general formulations in$^{44}$). We search for the eigenvalues $\mu_k^{(p)}$ and eigenfunctions $V_k^{(p)}$ of the generalized Steklov eigenvalue problem$^{31}$

$$ (p - \Delta) V_k^{(p)} = 0 \quad \text{in } \Omega, \quad \partial_n V_k^{(p)} = \mu_k^{(p)} V_k^{(p)} \quad \text{on } \partial\Omega, \quad \lim_{|x| \to \infty} V_k^{(p)} = 0, \qquad (2) $$

where $p > 0$ is a fixed parameter, $\Delta$ is the Laplace operator, and $\partial_n$ is the normal derivative oriented outwards the domain $\Omega$. The positive eigenvalues $\mu_k^{(p)}$ are ordered to form an increasing sequence that accumulates to infinity:

$$ 0 \le \mu_0^{(p)} \le \mu_1^{(p)} \le \dots \nearrow +\infty. \qquad (3) $$

We note that the restriction of $V_k^{(p)}$ onto the boundary $\partial\Omega$, denoted as $v_k^{(p)} = V_k^{(p)}|_{\partial\Omega}$, is an eigenfunction of the associated Dirichlet-to-Neumann operator $M_p$ that maps a given function $f$ on $\partial\Omega$ onto another function $g = (\partial_n u)|_{\partial\Omega}$, where $u$ is the solution of the boundary value problem

$$ (p - \Delta)u = 0 \quad \text{in } \Omega, \quad u|_{\partial\Omega} = f, \quad \lim_{|x| \to \infty} u = 0. \qquad (4) $$

As for bounded domains, the eigenfunctions $v_k^{(p)}$ form a complete orthonormal basis of the space $L^2(\partial\Omega)$ of square-integrable functions on the boundary $\partial\Omega^{31}$. In particular, their normalization also fixes the normalization of the Steklov eigenfunctions:

$$ \int_{\partial\Omega} |V_k^{(p)}|^2 = 1. \qquad (5) $$

In the conventional setting of bounded domains, solutions of Eq. (2) are searched in the space $H^1(\Omega)$ of square-integrable functions with square-integrable derivatives. When $p > 0$, solutions $V_k^{(p)}$ of Eq. (2) exhibit an exponential decay as $|x| \to \infty$ (see, e.g., upper bounds in Appendix A) and thus belong to the space $H^1(\Omega)$ (see$^{31}$ for further mathematical details). In turn, the analysis is more subtle for $p = 0$. A rigorous formulation of the exterior Steklov problem in two dimensions was provided in$^{45}$. In higher dimensions, one can either employ functional spaces larger than $H^1(\Omega)^{33,34,42}$, or consider the interior Steklov problem in a sequence of bounded domains approaching $\Omega^{43}$. In the following, we assume that $\mu_k^{(p)}$ and $V_k^{(p)}$ converge to their limits $\mu_k^{(0)}$ and $V_k^{(0)}$ that are *declared* to be eigenvalues and eigenfunctions of the exterior Steklov problem at $p = 0$. Mathematical proofs of this convergence and the equivalence with two other formulations of the exterior Steklov problem at $p = 0$ are beyond the scope of this paper and will be presented in a separate work$^{44}$.

In this paper, we aim at studying the limiting behavior of the Steklov eigenvalues and eigenfunctions as $p \to 0$ to highlight differences and similarities between interior and exterior Steklov problems. In particular, we address two questions:

(i) For a bounded domain, the principal eigenfunction $v_0^{(0)}$ corresponding to the smallest eigenvalue $\mu_0^{(0)} = 0$ is a constant. Does it hold for exterior domains?

(ii) For a bounded domain, the eigenvalues $\mu_k^{(p)}$ exhibit a linear asymptotic behavior as $p \to 0$,

$$ \mu_k^{(p)} = \mu_k^{(0)} + b_k p + o(p), \qquad (6) $$
where

$$b_k = \int_{\Omega} |V_k^{(0)}|^2 = \lim_{p \to 0} \partial_p \mu_k^{(p)} \qquad (7)$$

(see<sup>46</sup> and Appendix B). In particular, the smallest eigenvalue $\mu_0^{(p)}$ vanishes as $\mu_0^{(p)} \approx p|\Omega|/|\partial\Omega|$, where $|\Omega|$ and $|\partial\Omega|$ are the volume of $\Omega$ and the surface area of $\partial\Omega$, respectively, and $f \approx g$ means that the ratio $f/g$ tends to 1 in the considered limit. What is the small-$p$ asymptotic behavior of $\mu_k^{(p)}$ for exterior domains (for which $|\Omega| = \infty$)?

### III. EXAMPLE OF A BALL

To gain some intuition, it is instructive to start by considering the exterior of a ball $B_L \subset \mathbb{R}^d$ of radius $L$:

$$E_L = \{x \in \mathbb{R}^d : |x| > L\} = \mathbb{R}^d \setminus \overline{B_L}, \qquad (8)$$

where overline denotes the closure of $B_L$. The rotational symmetry of this domain and the consequent separation of variables in the modified Helmholtz equation written in spherical coordinates $(r, \theta, \dots)$ allow one to find the Steklov eigenvalues and eigenfunctions explicitly:

$$\mu_n^{(p,L)} = -\alpha \frac{k_{n,d}^{\prime}(\alpha L)}{k_{n,d}(\alpha L)}, \quad V_{n,m}^{(p,L)}(x) = \psi_{n,m}(\theta, \dots) g_n^{(p,L)}(r) \quad (n = 0, 1, 2, \dots), \qquad (9)$$

where $\alpha = \sqrt{p}$, $r = |x|$ is the radial coordinate, prime denotes the derivative with respect to the argument, $k_{n,d}(z) = z^{1-d/2}K_{n-1+d/2}(z)$ is an extension of the modified Bessel function $K_n(z)$ of the second kind, $\{\psi_{n,m}\}$ are normalized spherical harmonics on the sphere $\partial B_L$ in $\mathbb{R}^d$ (such that $\int_{\partial B_L} |\psi_{n,m}|^2 = 1$), and

$$g_n^{(p,L)}(r) = \frac{k_{n,d}(\alpha r)}{k_{n,d}(\alpha L)} \qquad (10)$$

are the radial functions. In contrast to the conventional enumeration of Steklov eigenfunctions in Eq. (2) by a single index $k$, we employ here a multi-index $n, m$ inherited from enumeration of spherical harmonics. For instance, the usual spherical harmonics $\psi_{n,m}(\theta, \phi)$ in three dimensions are enumerated by $n = 0, 1, 2, \dots$ and $m = -n, -n+1, \dots, n$. In two dimensions, one deals with Fourier harmonics $e^{in\theta}$ and $e^{-in\theta}$, which can alternatively be enumerated by $n = 0, 1, 2, \dots$ and $m = \pm 1$ distinguishing the sign. In higher dimensions, $m$ is a multi-index. As the eigenvalues $\mu_n^{(p,L)}$ do not depend on $m$, they are generally degenerate (e.g., $\mu_n^{(p,L)}$ is $(2n+1)$ times degenerate in three dimensions).

Setting $r = L$, one sees that spherical harmonics $\psi_{n,m}$ are the eigenfunctions of the Dirichlet-to-Neumann operator $\mathcal{M}_p^L$ on $\partial B_L$, associated with the eigenvalues $\mu_n^{(p,L)}$. In particular, these eigenfunctions do not depend on $p$; moreover, the principal eigenfunction, associated to the smallest eigenvalue $\mu_0^{(p,L)}$, is constant: $v_{0,0}^{(p,L)} = \psi_{0,0} = 1/\sqrt{|\partial B_L|}$.

One can easily check that the functions $k_{n,d}(z)$ satisfy the second-order differential equation

$$k_{n,d}''(z) + \frac{d-1}{z} k_{n,d}'(z) - \left(1 + \frac{n(n+d-2)}{z^2}\right) k_{n,d}(z) = 0, \qquad (11)$$

and the recurrence relation

$$k_{n,d}'(z) = -k_{n-1,d}(z) - \frac{n+d-2}{z}k_{n,d}(z), \qquad (12)$$

that yields

$$\mu_n^{(p,L)} = \frac{n+d-2}{L} + \alpha \frac{K_{n-2+d/2}(\alpha L)}{K_{n-1+d/2}(\alpha L)} \geq 0. \qquad (13)$$

Since $K_v(z)$ is a monotonously increasing function of $v$ for $z > 0$ and $v > 0$ (this property follows from an integral representation of $K_v(z)$, see 9.6.24 from<sup>47</sup>), we conclude that

$$0 \leq \mu_n^{(p,L)} - \mu_n^{(0,L)} \leq \sqrt{p} \qquad (14)$$
for any $n$ and $L$ in dimensions $d \ge 3$.

The exponential decay of the radial functions $g_n^{(p,L)}(r)$ as $r \to \infty$, inherited from the asymptotic behavior of $K_n(z)$ at large $z$, ensures the convergence of the $L^2(E_L)$-norm of all Steklov eigenfunctions $V_n^{(p,L)}$ for any $p > 0$, as expected$^{31}$. In Appendix C, we compute explicitly this norm and show that

$$ \int_{E_L} |V_{n,m}^{(p,L)}|^2 = \partial_p \mu_n^{(p,L)}. \qquad (15) $$

In the limit $p \to 0$, the asymptotic behavior of $K_n(z)$ implies

$$ \mu_n^{(0,L)} = \frac{n+d-2}{L}, \quad V_{n,m}^{(0,L)} = \psi_{n,m} (L/r)^{n+d-2} \quad (n=0,1,2,\dots). \qquad (16) $$

The small-$p$ asymptotic behavior of the norm $||V_{n,m}^{(p,L)}||_{L^2(E_L)}$ and of the associated eigenvalue $\mu_n^{(p,L)}$ turns out to be tightly related: if the $L^2(E_L)$-norm of $V_{n,m}^{(0,L)}$ is finite, one retrieves the asymptotic relation (6), derived for bounded domains; in turn, if the norm is infinite, a non-analytic term emerges in the leading order in $p$, as detailed below.

(i) In two dimensions, the limit of the principal eigenvalue is zero, $\mu_0^{(0,L)} = 0$, whereas the associated Steklov eigenfunction $V_{0,0}^{(p,L)}$ approaches a constant, $V_{0,0}^{(0,L)} = 1/\sqrt{2\pi L}$, which is not square integrable in the exterior of the disk. Similarly, the $L^2(E_L)$-norm of $V_{1,m}^{(0,L)}$ is also infinite. Accordingly, Eq. (13) implies an non-analytic behavior of the associated eigenvalues in the leading order (see also$^{27}$):

$$ \mu_0^{(p,L)} \approx \frac{-1}{L(\ln(\sqrt{pL}/2) + \gamma)} + O(p), \qquad (17a) $$

$$ \mu_1^{(p,L)} \approx \frac{1}{L} - Lp(\gamma + \ln(\sqrt{pL}/2)), \qquad (17b) $$

with $\gamma \approx 0.5772$ being the Euler's constant. In turn, the other eigenvalues exhibit a linear dependence on $p$ in the leading order:

$$ \mu_n^{(p,L)} \approx \frac{n}{L} + b_n^L p + o(p) \quad (n \ge 2), \qquad (18) $$

where

$$ b_n^L = \int_{E_L} |V_{n,m}^{(0,L)}|^2 = \frac{L}{2n+d-4} \qquad (19) $$

(see Appendix C for the computation of this integral).

(ii) In three dimensions, Eq. (13) yields

$$ \mu_0^{(p,L)} = \frac{1}{L} + \sqrt{p}, \quad \mu_n^{(p,L)} \approx \frac{n+1}{L} + b_n^L p + o(p) \quad (n \ge 1), \qquad (20) $$

where the first relation for $\mu_0^{(p,L)}$ is exact and valid for any $p \ge 0$. In other words, the asymptotic relation (6) fails for the principal eigenvalue, in agreement with the infinite $L^2(E_L)$-norm of $V_{0,0}^{(0,L)} = 1/(\sqrt{4\pi} |x|)$.

(iii) In four dimensions, the integral in Eq. (19) is infinite for $n = 0$, whereas the smallest eigenvalue $\mu_0^{(p,L)}$ behaves as

$$ \mu_0^{(p,L)} \approx \frac{2}{L} - Lp(\gamma + \ln(\sqrt{pL}/2)) + O(p). \qquad (21) $$

In turn, the other coefficients $b_n^L$ are finite and one retrieves the linear asymptotic behavior (6) for other eigenvalues.

(iv) Finally, for $d > 4$, the integral in Eq. (19) is finite for all eigenfunctions and Eq. (6) holds for all eigenvalues $\mu_n^{(p,L)}$. For instance, one has $\mu_0^{(p,L)} = 3/L + pL/(1 + \sqrt{pL}) \approx 3/L + pL + O(p^{3/2})$ for $d = 5$, i.e., the leading term is linear in $p$, whereas the non-analytic term appears in the next order.

In the next section, we extend the above results to other exterior domains and uncover the asymptotic behavior of the eigenvalues as $p \to 0$.
#### IV. ASYMPTOTIC BEHAVIOR OF EIGENVALUES

Let $V_k^{(p)}$ and $\mu_k^{(p)}$ be the Steklov eigenfunctions and eigenvalues in the exterior domain $\Omega = \mathbb{R}^d \setminus \Omega_0$, where $\Omega_0$ is a compact set in $\mathbb{R}^d$. We aim at deriving the asymptotic behavior of $\mu_k^{(p)}$ as $p \to 0$. For this purpose, we split $\Omega$ into a bounded domain $\Omega_L = \Omega \cap B_L$ and the exterior $E_L$ of a ball $B_L$ of radius $L$ large enough such that $\Omega_0 \subset B_L$. We start by deducing three identities.

##### A. Three identities

First, we integrate the equation $(p - \Delta) V_k^{(p)} = 0$ over the domain $\Omega_L$ to get

$$0 = p \int_{\Omega_L} V_k^{(p)} - \int_{\partial \Omega \cup \partial B_L} \partial_n V_k^{(p)}, \qquad (22)$$

so that

$$p \int_{\Omega_L} V_k^{(p)} = \mu_k^{(p)} \int_{\partial \Omega} v_k^{(p)} + \int_{\partial B_L} \partial_n V_k^{(p)}, \qquad (23)$$

where we used the Green's formula. As $V_k^{(p)}$ is an analytic function away from the boundary $\partial \Omega$, its normal derivative on $\partial B_L$ should be continuous:

$$(\partial_n V_k^{(p)})|_{\partial B_L^-} = -(\partial_n V_k^{(p)})|_{\partial B_L^+}, \qquad (24)$$

where $\partial B_L^- = \partial \Omega_L$ and $\partial B_L^+ = \partial E_L$ denote the spherical boundary $\partial B_L$ approached from inside and from outside, respectively, and the minus sign emerges from the opposite directions of the normal derivative on two sides of that boundary. The right-hand side of the above relation can then be expressed with the help of an auxiliary Dirichlet-to-Neumann operator $\mathcal{M}_p^L$ in the exterior of the ball $B_L$, implying

$$\partial_n V_k^{(p)} = - \mathcal{M}_p^L V_k^{(p)} \quad \text{on } \partial B_L \qquad (25)$$

(here we dropped the superscript minus for brevity). Integrating this expression over $\partial B_L$ and interpreting the integral as the scalar product in $L^2(\partial B_L)$, one gets

$$\int_{\partial B_L} \partial_n V_k^{(p)} = -(\mathcal{M}_p^L V_k^{(p)}, 1)_{L^2(\partial B_L)} = -(V_k^{(p)}, \mathcal{M}_p^L 1)_{L^2(\partial B_L)} = -\mu_0^{(p,L)} (V_k^{(p)}, 1)_{L^2(\partial B_L)} = -\mu_0^{(p,L)} \int_{\partial B_L} V_k^{(p)}, \qquad (26)$$

where $\mu_0^{(p,L)}$ is the principal eigenvalue of $\mathcal{M}_p^L$, associated to the constant eigenfunction (on the sphere $\partial B_L$). Substituting this expression into Eq. (23), we deduce the first identity

$$\mu_k^{(p)} \int_{\partial \Omega} v_k^{(p)} = p \int_{\Omega_L} V_k^{(p)} + \mu_0^{(p,L)} \int_{\partial B_L} V_k^{(p)}. \qquad (27)$$

In the limit $p \to 0$, one gets

$$\mu_k^{(0)} \int_{\partial \Omega} v_k^{(0)} = \frac{d-2}{L} \int_{\partial B_L} V_k^{(0)}. \qquad (28)$$

Second, we multiply $(p - \Delta) V_k^{(p)} = 0$ by $[V_k^{(p)}]^*$ and integrate over $\Omega_L$ to get

$$0 = \int_{\Omega_L} \left( p |V_k^{(p)}|^2 + |\nabla V_k^{(p)}|^2 \right) - \mu_k^{(p)} \underbrace{\int_{\partial \Omega} |v_k^{(p)}|^2}_{=1} - \int_{\partial B_L} [V_k^{(p)}]^* \partial_n V_k^{(p)}.$$
Using Eq. (25) and the eigenbasis of the operator $\mathcal{M}_p^L$, one obtains the following representation for the eigenvalue:

$$ \mu_k^{(p)} = \int_{\Omega_L} \left( p |V_k^{(p)}|^2 + |\nabla V_k^{(p)}|^2 \right) + \sum_{n,m} \mu_n^{(p,L)} |(V_k^{(p)}, \psi_{n,m})_{L^2(\partial B_L)}|^2. \quad (29) $$

In particular, one has in the limit $p \to 0$:

$$ \mu_k^{(0)} = \int_{\Omega_L} |\nabla V_k^{(0)}|^2 + \sum_{n,m} \mu_n^{(0,L)} |(V_k^{(0)}, \psi_{n,m})_{L^2(\partial B_L)}|^2. \quad (30) $$

This identity implies that all eigenvalues $\mu_k^{(0)}$ are strictly positive for $d \ge 3$:

$$ \mu_k^{(0)} > 0 \quad (k = 0, 1, 2, \dots). \quad (31) $$

Indeed, as $\mu_n^{(0,L)} = (n+d-2)/L > 0$ for all $n$ and $d \ge 3$, an eigenvalue $\mu_k^{(0)}$ could be zero only if all terms in Eq. (30) are strictly zero, i.e., if $V_k^{(0)}|_{\partial B_L}$ was orthogonal to all spherical harmonics $\psi_{n,m}$ that would contradict the completeness of $\{\psi_{n,m}\}$ in $L^2(\partial B_L)$. In other words, the condition $\mu_k^{(0)} = 0$ would require that $V_k^{(0)}|_{\partial B_L} \equiv 0$ for any $L$ large enough, that is impossible. Note that this argument is not applicable in two dimensions because $\mu_0^{(0,L)} = 0$, see Sec. IV B below.

Third, multiplying $(p-\Delta)V_k^{(p)} = 0$ by $[V_k^{(p')}]^*$, multiplying $(p'-\Delta)[V_k^{(p')}]^* = 0$ by $V_k^{(p)}$, subtracting these equations and integrating over $\Omega_L$, we have

$$ \begin{aligned} 0 &= (p-p') \int_{\Omega_L} V_k^{(p)} [V_k^{(p')}]^* - \int_{\Omega_L} ([V_k^{(p')}]^* \Delta V_k^{(p)} - V_k^{(p)} \Delta [V_k^{(p')}]^*) \\ &= (p-p') \int_{\Omega_L} V_k^{(p)} [V_k^{(p')}]^* - \int_{\partial \Omega \cup \partial B_L} ([V_k^{(p')}]^* \partial_n V_k^{(p)} - V_k^{(p)} \partial_n [V_k^{(p')}]^*) \\ &= (p-p') \int_{\Omega_L} V_k^{(p)} [V_k^{(p')}]^* - [\mu_k^{(p)} - \mu_k^{(p')}] \int_{\partial \Omega} v_k^{(p)} [v_k^{(p')}]^* + \int_{\partial B_L} \left( [V_k^{(p')}]^* \mathcal{M}_p^L V_k^{(p)} - V_k^{(p)} \mathcal{M}_{p'}^L [V_k^{(p')}]^* \right). \end{aligned} $$

We get thus the third identity:

$$ [\mu_k^{(p)} - \mu_k^{(p')}] \int_{\partial \Omega} v_k^{(p)} [v_k^{(p')}]^* = (p-p') \int_{\Omega_L} V_k^{(p)} [V_k^{(p')}]^* + W_k^{(p,p',L)}, \quad (32) $$

where

$$ \begin{aligned} W_k^{(p,p',L)} &= \int_{\partial B_L} \left( [V_k^{(p')}]^* \mathcal{M}_p^L V_k^{(p)} - V_k^{(p)} \mathcal{M}_{p'}^L [V_k^{(p')}]^* \right) = (V_k^{(p')}, (\mathcal{M}_p^L - \mathcal{M}_{p'}^L) V_k^{(p)})_{L^2(\partial B_L)} \\ &= \sum_{n,m} (V_k^{(p')}, \psi_{n,m})_{L^2(\partial B_L)} [\mu_n^{(p,L)} - \mu_n^{(p',L)}] (\psi_{n,m}, V_k^{(p)})_{L^2(\partial B_L)}. \end{aligned} \quad (33) $$

Here we used that the Dirichlet-to-Neumann operator $\mathcal{M}_{p'}^L$ is self-adjoint, while its eigenfunctions $\{\psi_{n,m}\}$ do not depend on $p$ and $p'$ (see Sec. III). Note that the scalar products in Eq. (33) can be expressed by using the identity (D3) from Appendix D.

The third identity has two important consequences. (i) Setting $p' = 0$, we get another identity that will be employed in the following analysis:

$$ [\mu_k^{(p)} - \mu_k^{(0)}] \int_{\partial \Omega} [v_k^{(0)}]^* v_k^{(p)} = p \int_{\Omega_L} [V_k^{(0)}]^* V_k^{(p)} + W_k^{(p,0,L)}. \quad (34) $$

(ii) Since the generalized exterior Steklov problem is well posed for $p > 0$, one can set $p' = p + \varepsilon$ and use a smooth dependence of $\mu_k^{(p)}$ and $V_k^{(p)}$ on $p$ and boundness of $\Omega_L$ to develop Eq. (32) in the limit $\varepsilon \to 0$

$$ \underbrace{\partial_p \mu_k^{(p)}}_{=1} \int_{\partial \Omega} |v_k^{(p)}|^2 = \int_{\Omega_L} |V_k^{(p)}|^2 + \sum_{n,m} |(V_k^{(p)}, \psi_{n,m})_{L^2(\partial B_L)}|^2 \partial_p \mu_n^{(p,L)}. \quad (35) $$
Substituting Eq. (15), one sees that the second term in the right-hand side is the integral of $|V_k^{(p)}|^2$ over the exterior $E_L$. Indeed, one has

$$ 
\begin{aligned}
\int_{E_L} |V_k^{(p)}|^2 &= \sum_{n,m,n'm'} (V_k^{(p)}, \psi_{n,m})_{L^2(\partial B_L)} (\psi_{n',m'}, V_k^{(p)})_{L^2(\partial B_L)} \int_{E_L} [V_{n,m}^{(p,L)}]^* V_{n',m'}^{(p,L)} \\
&= \sum_{n,m} |(V_k^{(p)}, \psi_{n,m})_{L^2(\partial B_L)}|^2 \underbrace{\int_{E_L} |V_{n,m}^{(p,L)}|^2}_{=\partial_p \mu_n^{(p,L)}}
\end{aligned}
$$

where the integral over angular coordinates canceled all the terms except if $n = n'$ and $m = m'$. As a consequence, we obtain

$$ \partial_p \mu_k^{(p)} = \int_{\Omega} |V_k^{(p)}|^2 \quad (p > 0). \qquad (36) $$

This relation, which was earlier derived for bounded domains$^{46}$, allows us to distinguish two scenarios in the limit $p \to 0$: if the squared $L^2(\Omega)$-norm of the limiting eigenfunction $V_k^{(0)}$ is finite, the associated eigenvalue $\mu_k^{(p)}$ approaches its limit $\mu_k^{(0)}$ linearly with $p$, as in Eq. (6) for bounded domains. In turn, if this norm is infinite, the derivative $\partial_p \mu_k^{(p)}$ is also infinite, suggesting an non-analytic behavior of $\mu_k^{(p)}$ near $p = 0$.

In the next subsection, we employ these identities to derive the asymptotic behavior of $\mu_k^{(p)}$ as $p \to 0$.

## B. Two-dimensional case

In two dimensions, the identity (28) implies for any $k$ that either $\mu_k^{(0)} = 0$, or the integral of $v_k^{(0)}$ over the boundary $\partial\Omega$ is zero. Christiansen and Datchev showed that $\mu_0^{(0)} = 0$, while $v_0^{(0)}$ is a constant function$^{45}$. According to Eq. (30), $\mu_0^{(0)} = 0$ implies that $V_0^{(0)}$ is also a constant function. Since such a constant function is unique (up to a multiplicative factor), the eigenvalue $\mu_0^{(0)}$ is simple, i.e., the next eigenvalue is strictly positive: $\mu_1^{(0)} > 0$. We conclude that

$$ \mu_0^{(0)} = 0, \qquad v_0^{(0)} = \frac{1}{\sqrt{|\partial\Omega|}}, \qquad (37) $$

$$ \mu_k^{(0)} > 0, \qquad \int_{\partial\Omega} v_k^{(0)} = 0 \quad (k > 0). \qquad (38) $$

The identity (34) yields in the leading order:

$$ \mu_k^{(p)} - \mu_k^{(0)} = \frac{\mu_0^{(p,L)}}{2\pi L} \left| \int_{\partial B_L} V_k^{(p)} \right|^2 + O(1/\ln^2(\sqrt{p})), \qquad (39) $$

where we wrote explicitly the leading-order term of $W_k^{(p,0,L)}$ from Eq. (33) with $n = 0$ (for which $\psi_{0,0} = 1/\sqrt{2\pi L}$), used the normalization of $v_k^{(0)}$ to replace $\int_{\partial\Omega} [v_k^{(0)}]^* v_k^{(p)}$ by 1 in the leading order, and ignored higher-order terms. Using Eq. (17a), we get in the leading order

$$ \mu_k^{(p)} - \mu_k^{(0)} = \frac{-a_k}{\ln(\sqrt{p})} + O(1/\ln^2(\sqrt{p})), \qquad (40) $$

where

$$ a_k = 2\pi|c_k|^2, \qquad c_k = \frac{1}{2\pi L} \int_{\partial B_L} V_k^{(0)} = \int_0^{2\pi} \frac{d\theta}{2\pi} V_k^{(0)}(L, \theta) \qquad (41) $$
(in the last equality, we used polar coordinates $(r, \theta)$ with the origin at the center of the disk $B_L$). It is important to stress that the coefficient $c_k$ does not depend on the choice of $L$. Indeed, the Steklov eigenfunction outside the disk of radius $L$ reads in polar coordinates $(r, \theta)$ as

$$V_k^{(0)}(r, \theta) = \sum_{n=-\infty}^{\infty} \left(\frac{L}{r}\right)^{|n|} e^{in\theta} \int_{0}^{2\pi} \frac{d\theta'}{2\pi} e^{-in\theta'} V_k^{(0)}(L, \theta') \quad (r \ge L). \quad (42)$$

One sees that the term with $n=0$ is precisely the coefficient $c_k$, which therefore determines a constant term in the asymptotic behavior of $V_k^{(0)}$ at large $|x|$:

$$V_k^{(0)}(x) \approx c_k + O(1/|x|) \quad |x| \to \infty. \quad (43)$$

When $c_k \ne 0$, the eigenvalue $\mu_k^{(p)}$ exhibits the logarithmically slow approach (40) to its limit $\mu_k^{(0)}$, while the associated eigenfunction approaches $c_k$ at infinity.

For instance, this asymptotic behavior holds for the principal eigenvalue with $k=0$. Indeed, as $v_0^{(0)} = 1/\sqrt{|\partial\Omega|}$, the associated eigenfunction $V_0^{(0)}$ is also constant, implying $c_0 = 1/\sqrt{|\partial\Omega|}$, and one gets $a_k = 2\pi/|\partial\Omega|$ and thus

$$\mu_0^{(p)} \approx \frac{-2\pi}{|\partial\Omega|\ln(\sqrt{p})} + O(1/\ln^2(\sqrt{p})) \quad (p \to 0). \quad (44)$$

Note that next-order terms determine the correct lengthscale that is needed to make the argument of the logarithm in the denominator dimensionless. Christiansen and Datchev analyzed the small-$p$ asymptotic behavior for connected unbounded domains $\Omega = \mathbb{R}^2 \setminus \Omega_0$ with a smooth boundary $\partial\Omega$ and identified this lengthscale as the logarithmic capacity $R_c$ of the compact set $\Omega_0$. In fact, they proved$^{45}$

$$\mu_0^{(p)} \approx \frac{-2\pi}{|\partial\Omega|(\ln(R_c\sqrt{p}/2) + \gamma)} + O(1/\ln^2(\sqrt{p})) \quad (p \to 0). \quad (45)$$

Even though this rigorous result has the same leading-order behavior as Eq. (44), it turns out to be more accurate, as will be checked by numerical simulations in Sec. V. Moreover, it provides a proper length $R_c$ to get a dimensionless argument of the logarithm.

For $k > 0$, we apply the identity (27) in the leading order,

$$\mu_k^{(p)} \int_{\partial\Omega} v_k^{(p)} = \mu_0^{(p,L)} \int_{\partial B_L} V_k^{(p)} + O(p), \quad (46)$$

as well as Eq. (17a, 40), to get

$$\mu_k^{(0)} \int_{\partial\Omega} v_k^{(p)} = \frac{-1}{\ln(\sqrt{p})} \left( \frac{1}{L} \int_{\partial B_L} V_k^{(0)} \right) + \dots, \quad (47)$$

where $\dots$ represents the next-order term. Using Eqs. (38, 41), we have

$$\int_{\partial\Omega} (v_k^{(p)} - v_k^{(0)}) = - \frac{2\pi c_k}{\mu_k^{(0)} \ln(\sqrt{p})} + O(1/\ln^2(\sqrt{p})) \quad (k > 0). \quad (48)$$

As a consequence, if $c_k \ne 0$, the eigenfunction $v_k^{(p)}$ approaches its limit $v_k^{(0)}$ logarithmically slowly,

$$v_k^{(p)} = v_k^{(0)} + \frac{u_k}{\ln(\sqrt{p})} + O(1/\ln^2(\sqrt{p})) \quad (k > 0), \quad (49)$$

where the leading-order correction term $u_k$ determines the coefficients $a_k$ and $c_k$:

$$c_k = - \frac{\mu_k^{(0)}}{2\pi} \int_{\partial\Omega} u_k, \quad a_k = \frac{[\mu_k^{(0)}]^2}{2\pi} \left| \int_{\partial\Omega} u_k \right|^2 \quad (k > 0). \quad (50)$$
This expression relates the asymptotic behavior of $V_k^{(0)}$ at infinity to the asymptotic behavior of $v_k^{(p)}$ on the boundary $\partial\Omega$.

What does happen in the case $c_k = 0$? The asymptotic behavior of $\mu_k^{(p)}$ is determined by the next-order term, which comes from the term $n = 1$ in Eq. (33). Indeed, the identity (34) yields in the leading order:

$$ \mu_k^{(p)} - \mu_k^{(0)} = \frac{\mu_1^{(p,L)}}{2\pi L} \left\{ \left| \int_{\partial B_L} V_k^{(0)} e^{i\theta} \right|^2 + \left| \int_{\partial B_L} V_k^{(0)} e^{-i\theta} \right|^2 \right\} + O(p), \quad (51) $$

where we wrote explicitly the Fourier harmonics $\psi_{1,m} = e^{im\theta}/\sqrt{2\pi L}$ in polar coordinates (with $m = \pm 1$). Substituting the leading-order term in the asymptotic relation (17b) for $\mu_1^{(p,L)}$, we find

$$ \mu_k^{(p)} - \mu_k^{(0)} = -p \ln(\sqrt{p}) d_k + O(p), \quad (52) $$

where

$$ d_k = \frac{L^2}{2\pi} \left\{ \left| \int_0^{2\pi} d\theta V_k^{(0)}(L, \theta) e^{i\theta} \right|^2 + \left| \int_0^{2\pi} d\theta V_k^{(0)}(L, \theta) e^{-i\theta} \right|^2 \right\}. \quad (53) $$

Using the identity (D3) derived in Appendix D, we can express the coefficient $d_k$ through the integral over the boundary $\partial\Omega$. In fact, Eq. (D3) reads in two dimensions for $n = 1$ and $m = \pm 1$:

$$ L \int_0^{2\pi} d\theta V_k^{(0)}(L, \theta) e^{im\theta} = \frac{1}{2} \int_{\partial\Omega} v_k^{(0)} \left( \mu_k^{(0)} re^{im\theta} - \partial_n(re^{im\theta}) \right), \quad (54) $$

where the polar coordinates $(r, \theta)$ were also used inside the disk $B_L$. We stress that both $r$ and $\theta$ in the right-hand side depend on the boundary point. As a consequence, the coefficient $d_k$ does not depend on $L$. While the coefficient $c_k$ determined the constant term in the asymptotic behavior of $V_k^{(0)}$ at infinity, the coefficient $d_k$ controls the $O(1/|x|)$ term.

Finally, if both $c_k$ and $d_k$ are zero, then $\mu_k^{(p)} - \mu_k^{(0)}$ in the identity (34) is controlled by $O(p)$ terms. In order to determine this asymptotic behavior, we note that both $O(1)$ and $O(1/|x|)$ terms in Eq. (42) are cancelled, i.e., $V_k^{(0)}$ behaves as $O(1/|x|^2)$ at infinity, so that its $L^2(\Omega)$ is finite. As a consequence, the limit of Eq. (36) as $p \to 0$ determines the derivative $\partial_p \mu_k^{(p)}$ at $p = 0$, and thus the asymptotic relation (6) holds, with $b_k$ given by Eq. (7), as for bounded domains.

The above findings can be summarized as follows:

(i) If $c_k \neq 0$, then $V_k^{(0)}(x) \simeq c_k$ as $|x| \to \infty$, and the asymptotic relation (40) holds;

(ii) If $c_k = 0$ and $d_k \neq 0$, then $V_k^{(0)} \simeq O(1/|x|)$ as $|x| \to \infty$, the $L^2(\Omega)$-norm of $V_k^{(0)}$ is infinite, and the asymptotic relation (52) holds;

(iii) If $c_k = d_k = 0$, then $V_k^{(0)} \simeq O(1/|x|^\nu)$ as $|x| \to \infty$ with $\nu \ge 2$, the $L^2(\Omega)$-norm of $V_k^{(0)}$ is finite, and the asymptotic relation (6) holds.

### C. Three-dimensional case

In three dimensions, one can substitute the asymptotic relations (20) into Eq. (33) to get in the leading order:

$$ W_k^{(p,0,L)} = \sqrt{p} |(V_k^{(0)}, \psi_{0,0})_{L^2(\partial B_L)}|^2 + O(p), \quad (55) $$

where $\psi_{0,0} = 1/\sqrt{4\pi L^2}$ is the constant eigenfunction of $\mathcal{M}_p^L$ on $\partial B_L$. As a consequence, the identity (34) implies in the leading order:

$$ [\mu_k^{(p)} - \mu_k^{(0)}] \underbrace{\int_{\partial\Omega} |v_k^{(0)}|^2}_{=1} = \frac{\sqrt{p}}{4\pi L^2} \left| \int_{\partial B_L} V_k^{(0)} \right|^2 + O(p). \quad (56) $$

Using the identity (28), we find

$$ \mu_k^{(p)} = \mu_k^{(0)} + a_k \sqrt{p} + O(p), \quad (57) $$
with

$$a_k = \frac{|\mu_k^{(0)}|^2}{4\pi} \left| \int_{\partial\Omega} v_k^{(0)} \right|^2. \qquad (58)$$

One sees that the integral of $v_k^{(0)}$ over the boundary $\partial\Omega$ determines the small-$p$ asymptotic behavior of the associated eigenvalue $\mu_k^{(p)}$. Note that Cauchy-Schwarz inequality implies an upper bound for the coefficients $a_k$:

$$a_k \le \frac{|\mu_k^{(0)}|^2}{4\pi} |\partial\Omega|. \qquad (59)$$

When the integral of $v_k^{(0)}$ over $\partial\Omega$ is zero, i.e., $a_k = 0$, one needs to search for the next-order term, which is linear in $p$. Substituting Eq. (20) for $n > 0$ into Eq. (33), we get

$$W_k^{(p,0,L)} = p \sum_{n,m} |(V_k^{(0)}, \psi_{n,m})_{L^2(\partial B_L)}|^2 \int_{E_L} |V_{n,m}^{(0,L)}|^2 + o(p), \qquad (60)$$

where there is no term with $n = 0$ (that yielded the $\sqrt{p}$ contribution in the above analysis) because $(V_k^{(0)}, \psi_{0,0})_{L^2(\partial B_L)} = 0$. Using Eqs. (C10, C12), one can evaluate the sum in Eq. (60) as

$$W_k^{(p,0,L)} = p \int_{E_L} |V_k^{(0)}|^2 + o(p). \qquad (61)$$

Substituting Eq. (61) into the identity (34), we get in the leading order:

$$[\mu_k^{(p)} - \mu_k^{(0)}] = p \int_{\Omega_L} |V_k^{(0)}|^2 + p \int_{E_L} |V_k^{(0)}|^2 + o(p) = p b_k + o(p),$$

where $b_k$ is the squared $L^2(\Omega)$-norm of $V_k^{(0)}$ defined in Eq. (7). We therefore obtained the asymptotic behavior (6), which was earlier derived for bounded domains.

We conclude that each eigenvalue $\mu_k^{(p)}$ of the generalized exterior Steklov problem in three dimensions exhibits one of two asymptotic relations:

$$\left\{ \begin{array}{ll} \text{(i) If } \int_{\partial\Omega} v_k^{(0)} \neq 0 & \Rightarrow a_k > 0, b_k = +\infty, \quad \mu_k^{(p)} = \mu_k^{(0)} + a_k \sqrt{p} + O(p), \\ \text{(ii) If } \int_{\partial\Omega} v_k^{(0)} = 0 & \Rightarrow a_k = 0, b_k < +\infty, \quad \mu_k^{(p)} = \mu_k^{(0)} + b_k p + o(p), \end{array} \right. \qquad (62)$$

where $a_k$ and $b_k$ are determined by Eqs. (58, 7), respectively.

For the exterior of a ball, the principal eigenfunction $v_0^{(0)}$ is constant so that other eigenfunctions $v_k^{(0)}$ are orthogonal to it, implying $a_k = 0$ for any $k > 0$ and thus the asymptotic behavior (6) for the eigenvalues $\mu_k^{(p)}$, in agreement with the explicit result (20). For a general exterior domain in $\mathbb{R}^3$, however, there is no reason for $v_0^{(0)}$ to be constant so that both asymptotic relations (6) and (57) are possible. For instance, the Steklov eigenfunctions for the exterior of prolate and oblate spheroids were shown to be non-constant$^{48}$. In Sec. V B, we solve numerically the generalized exterior Steklov problem for different domains and illustrate the predicted asymptotic behavior of $\mu_k^{(p)}$, as well as the pivotal role of the integral of $v_k^{(0)}$ over the boundary $\partial\Omega$.

We complete this subsection by considering the asymptotic behavior of $V_k^{(0)}$ at infinity (the following arguments are valid for $d \ge 3$). Let us assume that

$$V_k^{(0)} \approx C_k |x|^{-\gamma_k} \quad (|x| \to \infty), \qquad (63)$$

with some constants $C_k$ and $\gamma_k$. Substituting this behavior into the identity (28), one gets at large enough $L$:

$$\mu_k^{(0)} \int_{\partial\Omega} v_k^{(0)} \approx (d-2)S_d C_k L^{d-2-\gamma_k}, \qquad (64)$$
where $S_d = 2\pi^{d/2}/\Gamma(d/2)$ is the surface area of the unit sphere in $\mathbb{R}^d$. If the integral of $v_k^{(0)}$ over $\partial\Omega$ is not zero, one has

$$ \gamma_k = 2-d, \quad C_k = \frac{\mu_k^{(0)}}{(d-2)S_d} \int_{\partial\Omega} v_k^{(0)}. \qquad (65) $$

In contrast, if the integral of $v_k^{(0)}$ over $\partial\Omega$ is zero, the relation (64) cannot hold for any finite $L$, and the initial assumption (63) of an isotropic asymptotic behavior of $V_k^{(0)}$ fails.

### D. Higher-dimensional cases

In four dimensions, we use the asymptotic behavior (21) of the principal eigenvalue $\mu_0^{(p,L)}$ to get in the leading order

$$ W_k^{(p,0,L)} \approx -Lp \ln(\sqrt{p}) \frac{1}{2\pi^2 L^3} \left| \int_{\partial B_L} v_k^{(0)} \right|^2 + O(p), \qquad (66) $$

where $2\pi^2 L^3$ is the surface area of the sphere $\partial B_L$ in $\mathbb{R}^4$. From the identity (27), one can express the integral of $V_k^{(0)}$ over $\partial B_L$ to get

$$ W_k^{(p,0,L)} \approx -p \ln(\sqrt{p}) \frac{[\mu_k^{(0)}]^2}{8\pi^2} \left| \int_{\partial\Omega} v_k^{(0)} \right|^2 + O(p), \qquad (67) $$

and thus

$$ \mu_k^{(p)} = \mu_k^{(0)} - p \ln(\sqrt{p}) \frac{[\mu_k^{(0)}]^2}{8\pi^2} \left| \int_{\partial\Omega} v_k^{(0)} \right|^2 + O(p). \qquad (68) $$

Note that the linear term $O(p)$ determines the proper length scale to be inserted in $\ln(\sqrt{p})$, to make its argument dimensionless. As earlier, if the integral of $v_k^{(0)}$ over $\partial\Omega$ vanishes, then the leading term is linear in $p$.

When $d > 4$, one has $\mu_n^{(p,L)} - \mu_n^{(0,L)} \approx b_n^L p + o(p)$ for all Steklov eigenvalues for the exterior of a ball, so that the asymptotic relation (61) holds for all $k$, implying Eq. (6) for all Steklov eigenvalues $\mu_k^{(p)}$.

## V. NUMERICAL RESULTS

In this section, we undertake a numerical study of the generalized exterior Steklov problem in two and three dimensions. To our knowledge, the exterior of a ball is the only exterior domain, for which the eigenvalues and eigenfunctions of the Steklov problem are known explicitly (see Sec. III). We aim therefore to clarify how the shape of the exterior domain can affect the behavior of the Steklov eigenvalues and eigenfunctions. For this purpose, we adapt a finite-element method (FEM) that was developed for interior problems⁴⁹ by employing the transparent boundary condition to reduce the generalized exterior Steklov problem to an equivalent interior one (see details and numerical validation in Appendix E).

### A. Two dimensions

#### Eigenvalues

First, we numerically verify the logarithmic decay (45) of the principal eigenvalue $\mu_0^{(p)}$ for smooth shapes and simple polygons, for which the logarithmic capacity $R_c$ is analytically known⁵⁰⁻⁵²: (i) $R_c = R$ for a disk of radius $R$, (ii) $R_c = (a+b)/2$ for an ellipse with semiaxes $a$ and $b$, (iii) $R_c = \frac{\Gamma(1/4)^2}{4\pi^{3/2}} R$ for a square of sidelength $R$, and (iv) $R_c = \frac{\Gamma(1/3)^3}{2\pi^2\sqrt{3}} R$ for an equilateral triangle of sidelength $R$. Figure 1 presents the small-$p$ asymptotic behavior of $\mu_0^{(p)}$ for these four domains. Astonishingly, the asymptotic relation (45) turns out to be remarkably accurate, even though the correction term is formally of the order of $O(1/\ln^2(\sqrt{p}))$ (e.g.,
<figure><img src="image_2.png" /><figcaption>FIG. 1. The small-$p$ asymptotic behavior of $\mu_0^{(p)}$ for the exterior of (a) a disk of radius $R=1$; (b) an ellipse with semiaxes $a=1$ and $b=0.5$; (c) a square of sidelength $R=2$; and (d) an equilateral triangle of sidelength $R=2$. Solid line shows the asymptotic relation (45), while circles, asterisks and diamonds present the numerical results obtained with different mesh sizes and truncation orders (see Appendix E for details), that we used to highlight the numerical accuracy.</figcaption></figure>

$1/\ln(\sqrt{10^{-4}}) \approx -0.22$, i.e., even at very small $p$, the correction term could still be comparable to the leading one). Moreover, the asymptotic relation (45), which was derived in⁴⁵ for exterior domains with smooth boundary, works equally well for two considered polygonal domains, namely the square and the equilateral triangle. These results suggest that the asymptotic relation (45) can be more general (i.e., not limited to smooth shapes), while the correction term may be of higher order. An extension of Eq. (45) to more general domains with polygonal or Lipschitz boundary presents an interesting perspective.

According to our results from Sec. IV B, the other eigenvalues $\mu_k^{(p)}$ can approach their limits $\mu_k^{(0)}$ as $O(1/\ln(\sqrt{p}))$, $O(p\ln(\sqrt{p}))$, or $O(p)$. We observed all three types of these asymptotic relations for the considered examples of an ellipse, a square, and an equilateral triangle. Figure 2(a) presents the difference $\mu_k^{(p)} - \mu_k^{(0)}$ and the asymptotic relation (40) with $k=4$ (ellipse), $k=8$ (square), and $k=6$ (equilateral triangle). The coefficients $a_k$ from Eq. (41) were obtained by a numerical integration of $V_k^{(0)}$ over the circle of radius $L$. In turn, for other considered indices, the coefficients $a_k$ were very close to zero, suggesting a faster approach of $\mu_k^{(p)}$ to their limits $\mu_k^{(0)}$. In order to check this behavior, we present on panels (b,c,d) of Fig. 2 the ratio $(\mu_k^{(p)} - \mu_k^{(0)})/p$ for three earlier considered domains and several $k$. For the ellipse, we observe that $\mu_1^{(p)}$, $\mu_2^{(p)}$, and $\mu_5^{(p)}$ exhibit non-analytical small-$p$ behavior $O(p\ln(\sqrt{p}))$ according to the asymptotic relation (52), whereas the eigenvalue $\mu_3^{(p)}$ approaches to its limit linearly with $p$, via Eq. (6), as for bounded domains. After estimating the coefficients $d_k$ in front of the term $p\ln(\sqrt{p})$ via a numerical integration in Eq. (53), we plotted the asymptotic relation (52). These straight lines exhibit the correct slope (given by $d_k$). In turn, they miss the offsets, which are determined by the next-order terms $O(p)$ that we did not derive in Sec. IV B. A similar behavior is observed for the square (since $\mu_1^{(p)} = \mu_2^{(p)}$ is twice degenerate, the case $k=2$ is not shown). Note that the eigenvalues $\mu_3^{(p)}$ and $\mu_4^{(p)}$ are distinct but the associated coefficients $b_k$, determined by the squared $L^2(\Omega)$-norms of their Steklov eigenfunctions $V_k^{(0)}$, turn out to be identical. The curve for $k=5$ exhibits a small linear increase due to an non-analytical leading term $O(p\ln(\sqrt{p}))$ with a small numerical prefactor $d_5 \approx 0.05$. For the exterior of the equilateral triangle, the eigenvalues with $k=1,3,7$ exhibit an non-analytical behavior in the leading term, whereas the eigenvalue with $k=5$ behaves linearly with $p$. Note that we do not show the results for $k=2$ and $k=4$ because $\mu_2^{(p)} = \mu_1^{(p)}$ and $\mu_4^{(p)} = \mu_3^{(p)}$ due to the symmetries of the equilateral triangle.

## Eigenfunctions

We inspect the behavior of the principal eigenfunction $v_0^{(p)}$ of the Dirichlet-to-Neumann operator $\mathcal{M}_p$ for different shapes and values of $p$. Figure 3 illustrates the validity of Eq. (37) for a wide range of shapes such as the exterior of an ellipse, a square, an equilateral triangle and a polygonal shape shown in Fig. 4(bottom row). In contrast to the exterior of a disk, $v_0^{(p)}$ depends on $p$ for these shapes. Moreover, as $p$ increases, the eigenfunction exhibits variations of larger and larger amplitudes. For the last domain, $v_0^{(p)}$ becomes localized near the vertex with the smallest angle $\pi/6$. Note that the shown principal eigenfunctions $v_0^{(p)}$ are positive for all the considered cases. While this property is known for the interior Steklov problem$^{53,54}$, its extension to exterior domains remains an open question.

Since $v_0^{(0)}$ is constant, the related Steklov eigenfunction $V_0^{(0)}$ is also constant over the whole domain $\Omega$ (not shown). In turn, as soon as $p > 0$, the principal Steklov eigenfunction $V_0^{(p)}(x)$ exhibits an exponential decay at infinity (as $|x| \to \infty$). Panels (a) and (d) of Fig. 4 illustrate this behavior for the exterior of a square and of a polygonal domain. Two other Steklov eigenfunctions
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) The small-$p$ asymptotic behavior of $\mu_k^{(p)} - \mu_k^{(0)}$ with $k=4$ (circles) for the exterior of an ellipse with semiaxes $a=1$ and $b=0.5$; $k=8$ (asterisks) for the exterior of a square of sidelength $R=2$; and $k=6$ (diamonds) for the exterior of an equilateral triangle of sidelength $R=2$. Solid lines indicate the asymptotic relation (40), with $a_k$ being estimated as 0.062 (ellipse), 0.090 (triangle) and 0.047 (square). (b,c,d) The small-$p$ asymptotic behavior of $(\mu_k^{(p)} - \mu_k^{(0)})/p$ for the exterior of (b) an ellipse with semiaxes $a=1$ and $b=0.5$, (c) a square of sidelength $R=2$, and (d) an equilateral triangle of sidelength $R=2$. Symbols indicate the numerical results for different $k$, solid lines present the asymptotic relation (52) when $d_k \neq 0$, and dashed lines show the asymptotic relation (6) when $d_k = 0$. We have $d_1 \approx 0.62$, $d_2 \approx 0.82$, $d_3 = 0$ ($b_3 \approx 0.36$), $d_5 \approx 0.02$ (for ellipse); $d_1 \approx 1.03$, $d_3 = d_4 = 0$ ($b_3 \approx 0.55$, $b_4 \approx 0.56$), $d_5 \approx 0.05$ (for square); and $d_1 \approx 0.63$, $d_3 \approx 0.09$, $d_5 = 0$ ($b_5 \approx 0.21$), $d_6 \approx 0.05$ (for triangle). We used the maximal meshsize $h_{max} = 0.02$, the truncation order $n_{max} = 30$, and the outer boundary of radius $L=2$ (see Appendix E for details).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3. The principal eigenfunction $v_0^{(p)}$ with $p=0, 0.1, 1, 10$ for the exterior of (a) an ellipse with semiaxes $a=1$ and $b=0.5$; (b) a square $(-1, 1)^2$ of sidelength $R=2$; (c) an equilateral triangle of sidelength $R=2$; and (d) the polygonal shape shown in Fig. 4(bottom row). The horizontal axis corresponds to a curvilinear coordinate along the boundary $\partial\Omega$ (the index of equidistant boundary points), which follows the polar angle from $-\pi$ to $\pi$. For panels (a,b,c), the starting index corresponds to the boundary point with the polar angle $-\pi$, i.e., the point $(-1,0)$ for the ellipse and the square, the point $(-0.67,0)$ for the equilateral triangle, and the point $(1,-1)$ for the polygonal shape. Note that four and three spikes on panels (b) and (c) correspond to the vertices of the square and of the triangle, respectively. We used the maximal meshsize $h_{max} = 0.02$, the truncation order $n_{max} = 30$, and the outer boundary of radius $L=2$ (see Appendix E for details).</figcaption></figure>

$V_1^{(0)}$ and $V_3^{(0)}$ are also illustrated in Fig. 4. Both $V_0^{(0.01)}$ and $V_1^{(0)}$ are localized in the angle $\pi/6$ of the polygonal domain, in analogy to the interior problem (see$^{49}$ and references therein).

Figure 5 shows the decay of Steklov eigenfunctions $V_k^{(0)}$ for the exterior of a square. For this purpose, we plot $|V_k^{(0)}|$ along the negative part of the vertical axis. One sees that the eigenfunction with $k=1$ decays as $1/|x|$ so that their $L^2(\Omega)$-norm is infinite, yielding a non-analytic asymptotic behavior of the associated eigenvalues $\mu_k^{(p)}$, as shown in Fig. 2. In turn, the eigenfunction $V_3^{(0)}$ decays as $1/|x|^2$ so that its $L^2(\Omega)$-norm is finite, and the associated eigenvalue $\mu_3^{(p)}$ approaches linearly to its limit $\mu_3^{(0)}$. Finally, the eigenfunction $V_8^{(0)}$ does not decay but approaches a constant, as in Eq. (43), that agrees with the asymptotic behavior (40) of the eigenvalue $\mu_8^{(p)}$, as observed in Fig. 2(a).

## B. Three dimensions

A numerical analysis of the exterior Steklov problem in three dimensions can be carried out by the same finite-element method with the transparent boundary condition. However, a practical implementation of three-dimensional meshing and the consequent numerical diagonalization of large matrices can be rather time-consuming. For our illustrative purposes, it is sufficient to consider axisymmetric domains, for which the original three-dimensional problem can be treated as effectively two-dimensional (see
<figure><img src="image_5.png" /><figcaption>FIG. 4. Several Steklov eigenfunctions $V_k^{(p)}$ for the exterior of a square of sidelength $R = 2$ (top) and of a polygonal domain (bottom). We used the maximal meshsize $h_{max} = 0.02$, the truncation order $n_{max} = 30$, and the outer boundary of radius $L = 2$ (see Appendix E for details).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5. Decay of three Steklov eigenfunctions $|V_k^{(p)}|$ along a vertical line for the exterior of a square of sidelength $R = 2$, for three values of $p$. The associated eigenvalues are: $\mu_1^{(0)} \approx 0.76$, $\mu_3^{(0)} \approx 1.30$, $\mu_8^{(0)} \approx 3.37$; $\mu_1^{(0.01)} \approx 0.79$, $\mu_3^{(0.01)} \approx 1.31$, $\mu_8^{(0.01)} \approx 3.40$; and $\mu_1^{(1)} \approx 1.43$, $\mu_3^{(1)} \approx 1.70$, $\mu_8^{(1)} \approx 3.62$. A “spike” for the dashed curve corresponds to the sign change of the eigenfunction $V_8^{(0)}$. We used the maximal meshsize $h_{max} = 0.01$, the truncation order $n_{max} = 30$, and the outer boundary of radius $L = 2$ (see Appendix E for details).</figcaption></figure>

Appendix E for details). Moreover, we focus here exclusively on axisymmetric eigenfunctions, i.e., throughout this subsection, $V_k^{(p)}$ denotes the $(k+1)$-th axisymmetric Steklov eigenfunction, while $\mu_k^{(p)}$ is the associated eigenvalue. Other eigenvalues and eigenfunctions are not discussed (see$^{48}$ for a detailed study of various Steklov eigenfunctions with $p=0$ for prolate and oblate spheroids).

## Eigenvalues

Figure 6 illustrates the small-$p$ asymptotic behavior of first eigenvalues $\mu_k^{(p)}$ for the exterior of prolate/oblate spheroids, of a capped cylinder, and of a domain obtained rotating a given polygon around $z$ axis. For the case of prolate/oblate spheroids, the axisymmetric Steklov eigenfunctions $V_k^{(0)}$ were shown to be symmetric (resp., antisymmetric) with respect to the horizontal plane when $k$ is even (resp., odd)$^{48}$. As a consequence, the integral of $v_k^{(0)}$ over the boundary $\partial\Omega$ is zero for odd indices $k$,
<figure><img src="image_7.png" /><figcaption>FIG. 6. Small-$p$ asymptotic behavior of $\mu_k^{(p)} - \mu_k^{(0)}$ for the exterior of (a) a prolate spheroid with semiaxes $a = 0.5, b = 1$; (b) an oblate spheroid with semiaxes $a = 0.5, b = 1$; (c) a capped cylinder of radius 1 and height 2; and (d) a domain obtained by rotating the polygon shown in Fig. 7(bottom). Symbols show the numerical results (as indicated in the legend), whereas lines of the same color indicate the associated asymptotic behavior: solid lines for Eq. (57) and dashed lines for Eq. (6). A numerical computation of the coefficient $b_k$ given by Eq. (7) is described in Appendix C 1. We used the maximal meshsize $h_{max} = 0.02$, the truncation order $n_{max} = 30$, and the outer boundary of radius $L = 2$ (see Appendix E for details).</figcaption></figure>

implying $a_k = 0$ and the linear asymptotic relation (6). In turn, these integrals are nonzero for even indices $k$, yielding the asymptotic relation (57), in agreement with our numerical results shown on panels (a) and (b). Similarly, the symmetries of the capped cylinder imply that the integrals of $v_k^{(0)}$ over the boundary $\partial\Omega$ vanish for the eigenfunctions with $k = 1$ and $k = 3$. In turn, we obtain the asymptotic behavior (57) with $a_0 \approx 0.91$, $a_2 \approx 7.2 \cdot 10^{-4}$ and $a_4 \approx 0.07$ for the associated eigenvalues. Since the numerical value of $a_2$ is quite small, one needs to reach much smaller $p$ to ensure that $\sqrt{p}$ is indeed the leading term. In other words, the smallness of $a_2$ explains why the asymptotic relation (57) is not yet achieved at the considered values of $p$. Finally, for a general domain (without symmetries), there is no reason to get $a_k = 0$ for axisymmetric eigenfunctions, as illustrated on panel (d).

## Eigenfunctions

Top panels of Fig. 7 present first five axisymmetric Steklov eigenfunctions $V_k^{(0)}$ for the exterior of a capped cylinder of radius 1 and height 2. The shown projection on the $xz$ plane fully captures their behavior in the whole domain due to their axisymmetric symmetry. As for interior domains$^{53,54}$, the principal eigenfunction $V_0^{(0)}$ is positive, while the other eigenfunctions change their sign. One sees that the eigenfunctions $V_1^{(0)}$ and $V_3^{(0)}$ are antisymmetric with respect to the horizontal plane $xy$ (or, equivalently, to the horizontal axis for the shown projection) so that the integral of $v_k^{(0)}$ over the boundary $\partial\Omega$ is strictly zero, yielding $a_1 = a_3 = 0$, in agreement with the small-$p$ behavior shown in Fig. 6(c). In turn, the shown eigenfunctions $V_0^{(0)}, V_2^{(0)}$ and $V_4^{(0)}$ are symmetric with respect to the horizontal plane, suggesting that $a_k > 0$, in agreement with the previously reported values of $a_k$.

Bottom panels of Fig. 7 present first five axisymmetric Steklov eigenfunctions $V_k^{(0)}$ for the exterior of an axisymmetric domain obtained by rotating the shown polygon along the vertical axis. This polygon was chosen to have two obtuse angles $2\pi - \pi/3$ and $2\pi - \pi/6$ whose complementary angles $\pi/3$ and $\pi/6$ may accommodate localized Steklov eigenfunctions, e.g., $V_0^{(0)}$ on Fig. 7(f). Expectedly, a somewhat arbitrary shape of this polygon eliminates most symmetries (except for the axial symmetry imposed by construction) in the shown eigenfunctions. In particular, the integral of $v_k^{(0)}$ over $\partial\Omega$ is not zero for all indices $k \in \{0, 1, 2, 3, 4\}$, yielding $a_k > 0$ and thus the asymptotic behavior (57) shown in Fig. 6(d).

After examination of spatial variations of the Steklov functions near the boundary, we turn to their asymptotic behavior at large $|x|$. Figure 8 illustrates the expected power-law decay of Steklov eigenfunctions for $p = 0$ and the exponential decay (A10) of Steklov eigenfunctions for $p > 0$ for the exterior of a capped cylinder. Since $a_k > 0$ for $k = 0$ and $k = 2$, the associated eigenfunctions $|V_k^{(0)}|$ decreases as $1/|x|$, according to Eq. (63) with $\gamma_k = d - 2 = 1$. In turn, as $a_2 = 0$, the eigenfunction $V_2^{(0)}$ shows a faster power-law decay (here, with the exponent $-2$). When $p > 0$, all Steklov eigenfunctions exhibit an exponential decay, with an upper bound decreasing as $e^{-\sqrt{p}|x|}$ (see Appendix A). This behavior is confirmed on panels (b) and (c) of Fig. 8.
<figure><img src="image_8.png" /><figcaption>FIG. 7. First axisymmetric Steklov eigenfunctions $V_k^{(0)}$ for the exterior of (top) a capped cylinder of radius 1 and height 2; (bottom) an axisymmetric domain obtained by rotating the shown polygon along the vertical axis. The projection onto $xz$ plane, corresponding to the angular coordinate $\phi = 0$, is shown. We used the maximal meshsize $h_{max} = 0.02$, the truncation order $n_{max} = 30$, and the outer boundary of radius $L = 2$ (see Appendix E for details).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. Decay of three Steklov eigenfunctions $|V_k^{(p)}|$ along a vertical line for the exterior of a capped cylinder of radius 1 and height 2, for three values of $p$. Note that a spike for blue dash-dotted line ($k=2$) on panels (a) and (b) corresponds to a sign change of $V_k^{(p)}$. The associated eigenvalues are: $\mu_0^{(0)} \approx 0.78$, $\mu_1^{(0)} \approx 1.47$, $\mu_2^{(0)} \approx 1.96$; $\mu_0^{(0.01)} \approx 0.87$, $\mu_1^{(0.01)} \approx 1.48$, $\mu_2^{(0.01)} \approx 1.96$; $\mu_0^{(1)} \approx 1.70$, $\mu_1^{(1)} \approx 1.95$, $\mu_2^{(1)} \approx 2.27$. We used the maximal meshsize $h_{max} = 0.02$, the truncation order $n_{max} = 30$, and the outer boundary of radius $L = 2$ (see Appendix E for details).</figcaption></figure>

## VI. STATISTICS OF FIRST-PASSAGE TIMES

In this section, we describe an application of the small-$p$ asymptotic results to diffusion-controlled reactions and related first-passage time statistics. Indeed, as the modified Helmholtz equation $(p-\Delta)\tilde{u} = 0$ is related to the diffusion equation $\partial_t u = \Delta u$ via the Laplace transform, the small-$p$ behavior of $\tilde{u}(x, p)$ corresponds to the long-time behavior of $u(x, t)$. In particular, the encounter-based approach allowed one to express most characteristics of diffusion-reaction processes in the Laplace domain in terms of the Steklov eigenfunctions and eigenvalues$^{26–28,60–64}$.

To illustrate this point, we focus on the distribution of the first-crossing time $\mathcal{T}_l = \inf\{t > 0 : \ell_t > l\}$, i.e., the first time instance when the boundary local time $\ell_t$ crosses a fixed threshold $l$. The stochastic process $\ell_t$ naturally appears in the Skorokhod (or Langevin) stochastic equation of reflected Brownian motion in the presence of a reflecting boundary$^{55–59}$ and represents the
(rescaled) number of encounters between the diffusing particle and the boundary $\partial\Omega$ up to time $t$, and thus characterizes surface reactions$^{27}$. In turn, the random variable $\mathcal{T}_\ell$ plays the central role in defining the first-reaction times$^{28,60}$. The moment-generating function of the first-crossing time was expressed in$^{28}$ as

$$ \tilde{U}(\ell, p|x_0) = \langle e^{-p\mathcal{T}_\ell} \rangle = \int_0^\infty dt e^{-pt} U(\ell, t|x_0) = \sum_{k=0}^\infty v_k^{(p)}(x_0) e^{-\mu_k^{(p)}\ell} \int_{\partial\Omega} v_k^{(p)}, \quad (69) $$

where $x_0$ is the starting point, and $U(\ell, t|x_0)$ is the probability density function (PDF) of $\mathcal{T}_\ell$. For simplicity, we assume that the particle starts on the boundary $\partial\Omega$ (if $x_0 \notin \partial\Omega$, there is an additional stage of the first arrival onto $\partial\Omega$; as this first-passage process was thoroughly investigated in the past$^{8-11}$, we skip it here and set directly $x_0 \in \partial\Omega$).

Let us deduce the small-$p$ asymptotic behavior of Eq. (69) in three dimensions. In the leading order, one can replace $v_k^{(p)}$ by its limit $v_k^{(0)}$. As the resulting expression includes the integrals of $v_k^{(0)}$ over $\partial\Omega$, only the eigenmodes with nonzero integrals do contribute in the leading order. One can therefore use the asymptotic relation (57), with $a_k$ given by Eq. (58) so that

$$ \tilde{U}(\ell, p|x_0) \approx \sum_{k=0}^{\infty} v_k^{(0)}(x_0) e^{-(\mu_k^{(0)} + a_k \sqrt{p})\ell} \int_{\partial\Omega} v_k^{(0)} \quad (p \to 0), \quad (70) $$

from which the inverse Laplace transform yields

$$ U(\ell, t|x_0) \approx \sum_{k=0}^{\infty} \frac{a_k \ell}{\sqrt{4\pi Dt^3}} e^{-\mu_k^{(0)}\ell - a_k^2 \ell^2/(4Dt)} v_k^{(0)}(x_0) \int_{\partial\Omega} v_k^{(0)} \quad (t \to \infty), \quad (71) $$

where we included the diffusion coefficient $D$ (which was set to 1 in the remaining text). For instance, for a sphere of radius $L$, one has $\mu_0^{(0)} = 1/L$, $a_0 = 1$ and $v_0^{(0)} = 1/\sqrt{4\pi L^2}$, whereas the other eigenfunctions $v_k^{(0)}$ are orthogonal to $v_0^{(0)}$, yielding

$$ U(\ell, t|x_0) = \frac{\ell e^{-\ell/L - \ell^2/(4Dt)}}{\sqrt{4\pi Dt^3}} \quad (x_0 \in \partial\Omega). \quad (72) $$

Note that, in the case of the sphere, this expression turns out to be exact, i.e., it is valid for any $t$ (see$^{41}$).

The explicit form of the long-time asymptotic relation (71) highlights the relative roles of the eigenvalues $\mu_k^{(0)}$ and the coefficients $a_k$. In the common setting of a partially reactive boundary with a constant reactivity parameter $q$, the fixed threshold $\ell$ is replaced by a random threshold $\hat{\ell}$ with the exponential distribution, $\mathbb{P}\{\hat{\ell} > \ell\} = e^{-q\ell}$ (see$^{26,28}$). As a consequence, the probability density function $H_q(t|x_0)$ of the first-reaction time is obtained by averaging $U(\hat{\ell}, t|x_0)$ over all random realizations of $\hat{\ell}$:

$$ H_q(t|x_0) = \langle U(\hat{\ell}, t|x_0) \rangle = \int_0^\infty d\ell \underbrace{qe^{-q\ell}}_{\text{PDF of } \hat{\ell}} U(\ell, t|x_0). \quad (73) $$

Using the asymptotic relation (71), we get the long-time behavior (for $x_0 \in \partial\Omega$):

$$ H_q(t|x_0) \approx qD \sum_{k=0}^{\infty} \left\{ \frac{1}{\sqrt{\pi Dt}} - \frac{\mu_k^{(0)} + q}{a_k} \operatorname{erfcx} \left( \sqrt{Dt} \frac{\mu_k^{(0)} + q}{a_k} \right) \right\} \frac{v_k^{(0)}(x_0)}{a_k} \int_{\partial\Omega} v_k^{(0)}, \quad (74) $$

where $\operatorname{erfcx}(z) = e^{z^2} \operatorname{erfc}(z)$ is the scaled complementary error function. For the sphere of radius $L$, only the principal eigenmode contributes (with $a_0 = 1$ and $\mu_0^{(0)} = 1/L$), yielding

$$ H_q(t|x_0) = qD \left\{ \frac{1}{\sqrt{\pi Dt}} - (1/L + q) \operatorname{erfcx} \left( \sqrt{Dt} (1/L + q) \right) \right\} \quad (x_0 \in \partial\Omega), \quad (75) $$

and this asymptotic relation turns out to be exact (see$^{66}$ for details). To our knowledge, such a detailed representation (74) of the long-time asymptotic behavior of the PDF $H_q(t|x_0)$ was not earlier reported for general exterior domains in three dimensions.

At very long times, one can use the asymptotic behavior of $\operatorname{erfcx}(z)$ as $z \to \infty$ to get

$$ H_q(t|x_0) \approx \frac{q}{\sqrt{4\pi Dt^3}} \sum_{k=0}^{\infty} \frac{a_k}{(\mu_k^{(0)} + q)^2} v_k^{(0)}(x_0) \int_{\partial\Omega} v_k^{(0)}. \quad (76) $$
Expectedly, one retrieves the characteristic $t^{-3/2}$ heavy tail of the PDF in three dimensions, whereas the coefficients $a_k$ and eigenfunctions $v_k^{(0)}$ determine the prefactor.

Finally, the integral of $H_q(t|x_0)$ from $t$ to infinity determines the survival probability $S_q(t|x_0)$ of the particle in the presence of a partially reactive boundary $\partial\Omega$. Substituting Eq. (74) for large $t$, one gets the long-time asymptotic behavior of this survival probability (for $x_0 \in \partial\Omega$),

$$ S_q(t|x_0) \approx \sum_{k=0}^{\infty} \left\{ 1 - \frac{q}{\mu_k^{(0)} + q} \left[ 1 - \operatorname{erfcx} \left( \sqrt{Dt} \frac{\mu_k^{(0)} + q}{a_k} \right) \right] \right\} v_k^{(0)}(x_0) \int_{\partial\Omega} v_k^{(0)} \\ \approx 1 - q \sum_{k=0}^{\infty} \left[ 1 - \operatorname{erfcx} \left( \sqrt{Dt} \frac{\mu_k^{(0)} + q}{a_k} \right) \right] \frac{v_k^{(0)}(x_0)}{\mu_k^{(0)} + q} \int_{\partial\Omega} v_k^{(0)}, \quad (77) $$

where we used the completeness of the eigenfunctions $v_k^{(0)}$ to evaluate explicitly the first term (i.e., the constant 1). For the sphere, this asymptotic relation contains only the eigenmode with $k=0$ and turns out to be exact, reproducing the seminal result by Collins and Kimball$^{65}$. In the limit $t \to \infty$, the asymptotic relation (77) yields the exact form of the limiting survival probability (for $x_0 \in \partial\Omega$):

$$ S_q(\infty|x_0) = 1 - q \sum_{k=0}^{\infty} \frac{v_k^{(0)}(x_0)}{\mu_k^{(0)} + q} \int_{\partial\Omega} v_k^{(0)}. \quad (78) $$

This probability is nonzero due to the transient character of Brownian motion in three dimensions. In particular, $1 - S_q(\infty|x_0)$ is the probability of escape to infinity. Moreover, even though Eq. (77) was derived in the long-time limit, its right-hand side approaches the correct limit 1 as $t \to 0$. One can therefore expect that this relation might be a good approximation for the survival probability $S_q(t|x_0)$ over the whole range of times. A numerical validation of this approximation presents an interesting perspective.

Note that a similar analysis can be performed for the probability density function of the boundary local time$^{27,28}$. Moreover, one can also investigate the long-time asymptotic behavior of $U(\ell, t|x_0)$, $H_q(t|x_0)$ and $S_q(t|x_0)$ in two dimensions (see$^{41}$ for the related analysis in the case of the exterior of a disk).

## VII. CONCLUSION

In this paper, we studied the small-$p$ asymptotic behavior of the eigenvalues and eigenfunctions of the generalized exterior Steklov problem. The exponential decay of the Steklov eigenfunctions $V_k^{(p)}$ at infinity allowed us to explore the spectral properties for any $p > 0$ and then to consider the limit $p \to 0$. In particular, we derived the simple relation (36) between the squared $L^2(\Omega)$-norm of the eigenfunction $V_k^{(p)}$ and the derivative of the associated eigenvalue $\mu_k^{(p)}$ with respect to $p$. In the limit $p \to 0$, this relation distinguishes two situations: (i) when the $L^2(\Omega)$-norm of the limiting eigenfunction $V_k^{(0)}$ is finite, the associated eigenvalue $\mu_k^{(p)}$ approaches its limit $\mu_k^{(0)}$ linearly with $p$, as in the case of bounded domains; (ii) when the norm is infinite, an non-analytic behavior of $\mu_k^{(p)}$ is expected near $p = 0$. We managed to identify the leading non-analytic term for exterior domains in different space dimensions. In two dimensions, the principal eigenvalue $\mu_0^{(p)}$ vanishes logarithmically, as earlier proved rigorously in$^{45}$, while the associated eigenfunction becomes constant. This behavior was extended via Eq. (40) to other eigenvalues, whose eigenfunctions $V_k^{(0)}$ do not vanish at infinity (and thus the integrals $c_k$ defined by Eq. (41) were not zero). Numerical examples of such behavior were shown in Fig. 2(a). In turn, when $c_k = 0$ but the $L^2(\Omega)$-norm of $V_k^{(0)}$ is still infinite, the difference $\mu_k^{(p)} - \mu_k^{(0)}$ exhibits the non-analytic behavior $p \ln(\sqrt{p})$, which was observed for various exterior domains. Finally, if the $L^2(\Omega)$-norm of $V_k^{(0)}$ is finite, one retrieves the conventional linear relation (6), as for bounded domains. In three dimensions, the finiteness of the $L^2(\Omega)$-norm of $V_k^{(0)}$ is directly related to vanishing of the integral of the eigenfunction $v_k^{(0)}$ over the boundary $\partial\Omega$: if this integral is zero, the norm is finite, Eq. (36) holds, and $\mu_k^{(p)} - \mu_k^{(0)}$ is proportional to $p$ in the leading order. However, in the generic situation when the integral is not zero, the leading term becomes $\sqrt{p}$ in Eq. (57), with the coefficient $a_k$ given by Eq. (58). In the same vein, in four dimensions, if the integral of $v_k^{(0)}$ over $\partial\Omega$ is not zero, the non-analytic leading term $p \ln(\sqrt{p})$ emerges. In turn, all eigenfunctions $V_k^{(0)}$ have a finite $L^2(\Omega)$-norm in higher dimensions, implying the linearity of $\mu_k^{(p)} - \mu_k^{(0)}$ with $p$ in the leading order. These theoretical results were illustrated by several examples of exterior domains in two and three dimensions. For this purpose, we modified the earlier developed FEM by introducing the transparent boundary condition to deal with an equivalent Steklov problem in a bounded domain.
From the mathematical point of view, these asymptotic results clarify some properties of the exterior Steklov problem in the limit $p = 0$. For instance, this problem appears to be well posed in the conventional space $H^1(\Omega)$ in high dimensions $d \ge 5$. In turn, in dimensions $d = 2, 3, 4$, the $L^2(\Omega)$-norm of at least one limiting eigenfunction $V_k^{(0)}$ is expected to be infinite, thus requiring either a larger functional space$^{33,34,42}$, or a sequence of bounded domains approaching the unbounded one$^{43}$. A mathematical proof of the existence of the limit $p = 0$ and the inter-connections between various formulations of the exterior Steklov problem will be discussed in a separate work$^{44}$.

From the physical point of view, the small-$p$ asymptotic behavior of the eigenvalues and eigenfunctions of the generalized exterior Steklov problem has immediate consequences for diffusion-controlled reactions and related first-passage time statistics. We derived the long-time asymptotic behavior of the PDF $U(\ell, t|x_0)$ of the first-crossing time in three dimensions. Here the integrals of eigenfunctions $v_k^{(0)}$ over the boundary $\partial\Omega$ determine which Steklov eigenmodes do contribute to this asymptotic behavior, while the related coefficients $a_k$ control the exponential decay of these contributions. This PDF determines the statistics of various first-reaction times. In particular, we deduced the long-time asymptotic behavior of the survival probability and the related PDF of the first-reaction time in the common setting of a partially reactive boundary with constant reactivity. These asymptotic results emphasize the role of the coefficients $a_k$ and open a way to quantify the impact of the geometric shape of the boundary onto diffusion-controlled reactions in exterior domains.

## ACKNOWLEDGMENTS

The authors thank professors I. Polterovich and M. Levitin for fruitful discussions. D.S.G. acknowledges the Simons Foundation for supporting his sabbatical sojourn in 2024 at the CRM, University of Montréal, Canada. D.S.G. also acknowledges the Alexander von Humboldt Foundation for support within a Bessel Prize award.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## Appendix A: Upper bounds on Steklov eigenfunctions

In this Appendix, we inspect the decay of Steklov eigenfunctions at infinity. While this behavior is known, we employ some probabilistic arguments to get upper bounds.

As previously, we consider the exterior of a compact set $\Omega_0$, $\Omega = \mathbb{R}^d \setminus \Omega_0$, with a smooth boundary $\partial\Omega$. Each Steklov eigenfunction can be written as an extension of the corresponding eigenfunction $v_k^{(p)}$, defined on $\partial\Omega$:

$$V_k^{(p)}(x_0) = \int_{\partial\Omega} dx v_k^{(p)}(x) \tilde{j}_{\infty}(x, p|x_0), \quad (A1)$$

where $\tilde{j}_{\infty}(x, p|x_0) = \int_0^{\infty} dt e^{-pt} j_{\infty}(x, t|x_0)$ is the Laplace transform of the probability flux density $j_{\infty}(x, t|x_0) = -\partial_n G_{\infty}(x, t|x_0)$ (hereafter tilde denotes the Laplace transform). Here $G_{\infty}(x, t|x_0)$ is the heat kernel (or propagator), satisfying the diffusion equation, $\partial_t G_{\infty}(x, t|x_0) = \Delta G_{\infty}(x, t|x_0)$, with Dirichlet condition, $G_{\infty}(x, t|x_0)|_{\partial\Omega} = 0$, and the initial condition $G_{\infty}(x, 0|x_0) = \delta(x - x_0)$ with a Dirac distribution $\delta(x - x_0)$. Note that the diffusion coefficient is set to 1 for brevity.

As a consequence, one gets

$$|V_k^{(p)}(x_0)| \le \|v_k^{(p)}\|_{L^{\infty}(\partial\Omega)} \tilde{H}_{\infty}(p|x_0), \quad (A2)$$

where

$$\tilde{H}_{\infty}(p|x_0) = \int_{\partial\Omega} dx \tilde{j}_{\infty}(x, p|x_0) \quad (A3)$$

is the Laplace transform of the probability density $H_{\infty}(t|x_0)$ of the first-passage time $\tau$ to $\partial\Omega$, whereas $\|v_k^{(p)}\|_{L^{\infty}(\partial\Omega)}$ is the maximum of $|v_k^{(p)}|$ over $\partial\Omega$. Note that the probability density $H_{\infty}(t|x_0)$ is not normalized to 1 for $d \ge 3$ due to the transient character of Brownian motion.
Let $\Omega_0 \subset \Omega_0^*$ so that $\Omega^* = \mathbb{R}^d \setminus \Omega_0^* \subset \Omega$. For a starting point $x_0 \in \Omega^*$, we introduce the first-passage time $\tau^*$ to the boundary $\partial\Omega^*$. Since the trajectories of Brownian motion are continuous, each trajectory arriving to $\partial\Omega$ must cross $\partial\Omega^*$. As a consequence, one can formally write $\tau^* \le \tau$, which means that

$$S_{\infty}^{*}(t|x_0) = \mathbb{P}_{x_0}\{t < \tau^{*}\} \le \mathbb{P}_{x_0}\{t < \tau\} = S_{\infty}(t|x_0), \quad (A4)$$

where $S_{\infty}(t|x_0)$ and $S_{\infty}^{*}(t|x_0)$ are the survival probabilities in $\Omega$ and $\Omega^*$, respectively. Multiplying this inequality by $e^{-pt}$ and integrating over $t$ from 0 to infinity, one gets for the Laplace transforms:

$$\tilde{S}_{\infty}^{*}(p|x_0) \le \tilde{S}_{\infty}(p|x_0) \quad (p \ge 0). \quad (A5)$$

In turn, as $H_{\infty}(t|x_0) = -\partial_t S_{\infty}(t|x_0)$ and $H_{\infty}^{*}(t|x_0) = -\partial_t S_{\infty}^{*}(t|x_0)$, one has

$$\tilde{H}_{\infty}(p|x_0) = 1 - p\tilde{S}_{\infty}(p|x_0) \le 1 - p\tilde{S}_{\infty}^{*}(p|x_0) = \tilde{H}_{\infty}^{*}(p|x_0). \quad (A6)$$

We conclude that $\tilde{H}_{\infty}^{*}(p|x_0)$ plays the role of an upper bound for $\tilde{H}_{\infty}(p|x_0)$ and thus for $|V_k^{(p)}(x_0)|$:

$$|V_k^{(p)}(x_0)| \le \|v_k^{(p)}\|_{L^{\infty}(\partial\Omega)} \tilde{H}_{\infty}^{*}(p|x_0). \quad (A7)$$

Choosing $\Omega_0^* = \{x \in \mathbb{R}^d : |x| < L\}$ to be a ball of radius $L$ that englobes $\Omega_0$, one has

$$\tilde{H}_{\infty}^{*}(p|x_0) = (|x_0|/L)^{1-\frac{d}{2}} \frac{K_{\frac{d}{2}-1}(|x_0|\sqrt{p})}{K_{\frac{d}{2}-1}(L\sqrt{p})}, \quad (A8)$$

where $K_v(z)$ is the modified Bessel function of the second kind. For $p > 0$ and large $|x_0|$, one uses the asymptotic behavior of $K_v(z)$ to get

$$\tilde{H}_{\infty}^{*}(p|x_0) \approx (|x_0|/L)^{1-\frac{d}{2}} \sqrt{\frac{\pi}{2} K_{\frac{d}{2}-1}^2(L\sqrt{p})} e^{-|x_0|\sqrt{p}}. \quad (A9)$$

As a consequence, each Steklov eigenfunction $V_k^{(p)}(x_0)$ decays exponentially fast at infinity for any $p > 0$:

$$|V_k^{(p)}(x_0)| \le C_{k,p} e^{-|x_0|\sqrt{p}} \quad (|x_0| \gg L), \quad (A10)$$

with an non-universal constant $C_{k,p} > 0$. In particular, this property ensures that the $L^2(\Omega)$-norm of $V_k^{(p)}$ is finite for any $p > 0$. In turn, for $p = 0$, one gets

$$\tilde{H}_{\infty}^{*}(0|x_0) = \begin{cases} (|x_0|/L)^{2-d} & (d \ge 3), \\ 1 & (d = 2). \end{cases} \quad (A11)$$

This is the hitting probability of the ball of radius $L$. Due to the recurrent nature of Brownian motion in two dimensions, this probability is equal to 1; in turn, it is smaller than 1 for Brownian motion in higher dimensions due to the possibility of escape to infinity (transient motion). As a consequence, each Steklov eigenfunction $V_k^{(0)}(x_0)$ is bounded by an explicit function that decays via a power law as $|x_0| \to \infty$ in dimensions $d \ge 3$. In turn, the above upper bound is useless for $d = 2$.

## Appendix B: Small-$p$ asymptotic behavior for bounded domains

For a bounded domain, the small-$p$ asymptotic behavior of the eigenvalues were discussed earlier, see, e.g.,$^{27,67}$. In$^{49}$, we used probabilistic arguments to get the asymptotic relation (6). In this conventional setting, the squared $L^2(\Omega)$-norm of $V_k^{(0)}$, given by $b_k$, is always finite, ensuring the asymptotic behavior (6). In this Appendix, we briefly discuss a direct derivation of this result for a more general mixed Steklov problem for bounded domains.

We consider a bounded domain $\Omega \subset \mathbb{R}^d$ whose boundary $\partial\Omega$ is split into two parts, $\partial\Omega_S \cup \partial\Omega_D$, with Steklov and Dirichlet boundary conditions on $\partial\Omega_S$ and $\partial\Omega_D$, respectively. As in Sec. IV A, we multiply $(p - \Delta) V_k^{(p)} = 0$ by $[V_k^{(p')}]^*$, multiply $(p' - \Delta)[V_k^{(p')}]^* = 0$ by $V_k^{(p)}$, subtract these equations and integrate over $\Omega$ to get

$$\begin{aligned} 0 &= (p - p') \int_{\Omega} V_k^{(p)} [V_k^{(p')}]^* - \int_{\Omega} ([V_k^{(p')}]^* \Delta V_k^{(p)} - V_k^{(p)} \Delta [V_k^{(p')}]^*) \\ &= (p - p') \int_{\Omega} V_k^{(p)} [V_k^{(p')}]^* - \int_{\partial\Omega_S} ([V_k^{(p')}]^* \partial_n V_k^{(p)} - V_k^{(p)} \partial_n [V_k^{(p')}]^*) - \int_{\partial\Omega_D} ([V_k^{(p')}]^* \partial_n V_k^{(p)} - V_k^{(p)} \partial_n [V_k^{(p')}]^*). \end{aligned} \quad (B1)$$
The last term vanishes due to the Dirichlet boundary condition on $\partial\Omega_D$, yielding

$$ (p' - p) \int_{\Omega} V_k^{(p)} [V_k^{(p')}]^* = [\mu_k^{(p')} - \mu_k^{(p)}] \int_{\partial\Omega_S} v_k^{(p)} [v_k^{(p')}]^*. \quad (B2) $$

From a general theory of elliptic equations in bounded domains⁹¹, the eigenvalues exhibit analytical dependence on $p$ (for any $p > 0$), up to their suitable re-ordering to stay on the chosen analytic branch (here we exclude the case of degenerate eigenvalues that requires more subtle analysis, see also⁹²,⁹³). Setting $p' = p + \epsilon$ and taking the limit $\epsilon \to 0$, one gets

$$ \int_{\Omega} |V_k^{(p)}|^2 = \partial_p \mu_k^{(p)} \quad (p \ge 0). \quad (B3) $$

As the $L^2(\Omega)$-norm of $V_k^{(0)}$ is finite for a bounded domain, one retrieves Eq. (6), with the coefficients $b_k$ given by Eq. (7). Note that this derivation remains valid for mixed Steklov-Neumann, Steklov-Robin and Steklov-Robin-Dirichlet problems when the Dirichlet boundary condition on $\partial\Omega_D$ is replaced by Neumann or Robin boundary condition. Indeed, for all these boundary conditions, the last term in Eq. (B1) still vanishes.

## Appendix C: Explicit computation for the exterior of a ball

In this Appendix, we focus on the exterior of a ball of radius $L$, $E_L = \mathbb{R}^d \setminus \overline{B_L}$, and evaluate the integral

$$ Q_n^{(p,L)} = \int_{E_L} |V_{n,m}^{(p,L)}|^2 = \int_{E_L} |\psi_{n,m}|^2 [g_n^{(p,L)}(r)]^2 = \frac{1}{L^{d-1}} \int_L dr r^{d-1} [g_n^{(p,L)}(r)]^2, \quad (C1) $$

where we used the factorization (9) of a Steklov eigenfunction $V_{n,m}^{(p,L)}$ as the product of a spherical harmonic $\psi_{n,m}$ and a radial function $g_n^{(p,L)}$. Note that the integral over angular coordinates yielded the factor $1/L^{d-1}$ due to the normalization of spherical harmonics: $\int_{\partial B_L} |\psi_{n,m}|^2 = 1$. Using the explicit form (10) of the radial function, one has

$$ Q_n^{(p,L)} = \frac{\mathcal{I}_n(\alpha L)}{L^{d-1} [k_{n,d}(\alpha L)]^2 \alpha^d}, \quad (C2) $$

where $\alpha = \sqrt{p}$, and

$$ \mathcal{I}_n(z) = \int_z^\infty dx x^{d-1} [k_{n,d}(x)]^2. \quad (C3) $$

The last integral can be found in a standard way by using the differential equation (11) for $y(x) = k_{n,d}(x)$. Multiplying this equation by $y(x)$ and integrating from $z$ to $\infty$, one gets

$$ \int_z^\infty dx x^{d-1} \left( [y'(x)]^2 + [y(x)]^2 \left(1 + n(n+d-2)/x^2\right) \right) = -z^{d-1}y(z)y'(z). \quad (C4) $$

In turn, multiplying Eq. (11) by $xy'(x)$ and integrating from $z$ to $\infty$, one has

$$ \begin{aligned} & \int_z^\infty dx x^{d-1} \left( (d-2)[y'(x)]^2 + [y(x)]^2 \left(d + (d-2)n(n+d-2)/x^2\right) \right) \\ & = z^d[y'(z)]^2 - z^d[y(z)]^2 \left(1 + n(n+d-2)/z^2\right). \end{aligned} \quad (C5) $$

Multiplying Eq. (C4) by $d-2$ and subtracting from Eq. (C5), one finds

$$ \mathcal{I}_n(z) = \frac{z^{d-2}}{2} \left( (d-2)zy(z)y'(z) + z^2[y'(z)]^2 - y^2(z)[z^2 + n(n+d-2)] \right). \quad (C6) $$
Using the recurrence relation (12) we further simplify this expression as

$$ \mathcal{J}_n(z) = \frac{z^{d-1} k_{n,d}^2(z)}{2} \left( (2n+d-2) \frac{k_{n-1,d}(z)}{k_{n,d}(z)} + z \left[ \frac{k_{n-1,d}^2(z)}{k_{n,d}^2(z)} - 1 \right] \right), \quad (C7) $$

and thus

$$ Q_n^{(p,L)} = \frac{1}{2\alpha} \left( (2n+d-2) \frac{k_{n-1,d}(\alpha L)}{k_{n,d}(\alpha L)} + \alpha L \left[ \frac{k_{n-1,d}^2(\alpha L)}{k_{n,d}^2(\alpha L)} - 1 \right] \right). \quad (C8) $$

At the same time, one can easily check that the derivative of $\mu_n^{(p,L)}$ given by Eq. (9) with respect to $p$ yields the same expression. We conclude that

$$ Q_n^{(p,L)} = \partial_p \mu_n^{(p,L)}, \quad (C9) $$

i.e., the squared $L^2(E_L)$-norm of $V_{n,m}^{(p,L)}$ is determined by the derivative of the associated eigenvalue. We therefore retrieved the relation (7) derived earlier for bounded domains.

In the limit $p \to 0$, one uses the asymptotic behavior of the modified Bessel function to get

$$ k_{n,d}(z) \approx \Gamma(n-1+d/2) 2^{n-2+d/2} z^{-(n+d-2)} \quad (z \to 0) $$

so that $k_{n-1,d}(z)/k_{n,d}(z) \approx z/(2n+d-4)$ for $2n+d > 4$ and thus

$$ Q_n^{(0,L)} = \begin{cases} L/(2n+d-4) & (n > 2-d/2), \\ +\infty & (n \le 2-d/2). \end{cases} \quad (C10) $$

One can easily check that the divergence for $2n+d \le 4$ occurs at $n=0$ for dimensions $d=2,3,4$, and at $n=1$ for dimension $d=2$.

We use this result for evaluating the integral of a harmonic function $u$ in the exterior of the ball. Writing

$$ u = \sum_{n,m} V_{n,m}^{(0,L)} (u|_{\partial B_L}, \psi_{n,m})_{L^2(\partial B_L)}, \quad (C11) $$

one has

$$ \int_{E_L} |u|^2 = \sum_{n_1,m_1,n_2,m_2} (u|_{\partial B_L}, \psi_{n_1,m_1})_{L^2(\partial B_L)}^* (u|_{\partial B_L}, \psi_{n_2,m_2})_{L^2(\partial B_L)} \int_{E_L} [V_{n_1,m_1}^{(0,L)}]^* V_{n_2,m_2}^{(0,L)} \\ = \sum_{n,m} |(u|_{\partial B_L}, \psi_{n,m})_{L^2(\partial B_L)}|^2 \int_{E_L} |V_{n,m}^{(0,L)}|^2, $$

where we used the orthogonality of spherical harmonics $\{\psi_{n,m}\}$. As a consequence, we obtain

$$ \int_{E_L} |u|^2 = \sum_{n,m} Q_n^{(0,L)} |(u|_{\partial B_L}, \psi_{n,m})_{L^2(\partial B_L)}|^2. \quad (C12) $$

For instance, in three dimensions, if the integral of $u|_{\partial B_L}$ over the sphere $\partial B_L$ is zero, the divergent term $Q_0^{(0,L)}$ is excluded, and the integral (C12) is finite.

## 1. Numerical evaluation of the coefficient $b_k$

To compute the coefficients $b_k$ from Eq. (7), one can split the unbounded domain $\Omega$ into $\Omega_L = \Omega \cap B_L$ and the exterior of the ball $B_L$, $E_L = \mathbb{R}^d \setminus B_L$ so that

$$ b_k = \int_{\Omega_L} |V_k^{(0)}|^2 + \int_{E_L} |V_k^{(0)}|^2. \quad (C13) $$
The first term can be calculated directly by using the FEM solution inside the computational domain $\Omega_L$. As discussed earlier, if the integral of $v_k^{(0)}$ over $\partial\Omega$ is not zero, the coefficient $b_k$ is infinite. We therefore focus on the case when this integral is zero and restrict the following discussion to three dimensions. According to Eq. (28), this implies that $(V_k^{(0)}, 1)_{L^2(\partial B_L)} = 0$. As $\psi_{0,0} = 1/\sqrt{4\pi L^2}$ is constant, one can exclude the term with $n=0$ from the sum in Eq. (C12), yielding

$$ \int_{E_L} |V_k^{(0)}|^2 = \sum_{n>0} \sum_{m=-n}^{n} \left| (\psi_{n,m}, V_k^{(0)})_{L^2(\partial B_L)} \right|^2 \frac{L}{2n-1} \quad (C14) $$

in three dimensions.

If we focus on the axisymmetric Steklov eigenfunctions, which do not depend on the angle $\phi$, one can employ axisymmetric spherical harmonics given by (normalized) Legendre polynomials:

$$ \psi_{n,0}(\theta) = \sqrt{\frac{2n+1}{4\pi L^2}} P_n(\cos \theta). \quad (C15) $$

As a consequence, we have

$$ (\psi_{n,0}, V_k^{(0)})_{L^2(\partial B_L)} = L^2 \int_0^{2\pi} d\phi \int_0^{\pi} d\theta \sin\theta \psi_{n,0}(\theta) V_k^{(0)}(L, \theta), \quad (C16) $$

so that

$$ \int_{E_L} |V_k^{(0)}|^2 = 2\pi L^3 \sum_{n=1}^{\infty} \frac{2n+1}{2n-1} \left| \int_0^{\pi} d\theta \sin\theta P_n(\cos\theta) V_k^{(0)}(L, \theta) \right|^2. \quad (C17) $$

In practice, this integral can be approximated by a Riemann sum which is then evaluated over the nodes of the sphere $\partial B_L$, see Appendix E.

In two dimensions, the $L^2(\Omega)$-norm of $V_k^{(0)}$ is finite only if the first two terms with $n=0$ and $n=1$ are excluded (i.e., their coefficients are strictly zero). In this case, we get

$$ \begin{aligned} \int_{E_L} |V_k^{(0)}|^2 &= \sum_{n \ge 2} \frac{L}{2n-2} \sum_{m=\pm 1} \left| (\psi_{n,m}, V_k^{(0)})_{L^2(\partial B_L)} \right|^2 \\ &= \frac{L^2}{4\pi} \sum_{n \ge 2} \frac{1}{n-1} \left\{ \left| \int_0^{2\pi} d\theta e^{in\theta} V_k^{(0)}(L, \theta) \right|^2 + \left| \int_0^{2\pi} d\theta e^{-in\theta} V_k^{(0)}(L, \theta) \right|^2 \right\}, \end{aligned} \quad (C18) $$

where we used $\psi_{n,m}(\theta) = e^{imn\theta}/\sqrt{2\pi L}$ (alternatively, one can use $\cos(n\theta)$ and $\sin(n\theta)$).

## Appendix D: Fourth identity

Here we derive another identity, which can be used to evaluate and analyze the terms in Eq. (33). For this purpose, we consider the Steklov eigenvalues $\hat{\mu}_{n,m}^{(p,L)}$ and eigenfunctions $\hat{V}_{n,m}^{(p,L)}$ for the interior of a ball of radius $L$:

$$ \hat{\mu}_n^{(p,L)} = \alpha \frac{I_{n-2+d/2}(\alpha L)}{I_{n-1+d/2}(\alpha L)} - \frac{n+d-2}{L}, \quad \hat{V}_{n,m}^{(p,L)}(x) = \psi_{n,m}(\theta, \dots) \frac{i_{n,d}(\alpha r)}{i_{n,d}(\alpha L)} \quad (n=0,1,2,\dots), \quad (D1) $$

where $\alpha = \sqrt{p}$, $i_{n,d}(z) = z^{1-d/2}I_{n-1+d/2}(z)$ is an extension of the modified Bessel function $I_v(z)$ of the first kind, and we used the recurrence relation

$$ i'_{n,d}(z) = i_{n-1,d}(z) - \frac{n+d-2}{z} i_{n,d}(z). \quad (D2) $$

In the limit $p \to 0$, one gets $\hat{\mu}_n^{(0,L)} = n/L$ and $\hat{V}_{n,m}^{(0,L)}(x) = \psi_{n,m}(\theta, \dots)(r/L)^n$; in particular, the principal eigenvalue is zero, as for any interior problem.
Using the eigenfunctions $\hat{V}_{n,m}^{(p,L)}$, we can express the integrals of $V_k^{(p)}$ multiplied by $\psi_{n,m}$ on the sphere $\partial B_L$. Indeed, multiplying $(p-\Delta)V_k^{(p)} = 0$ by $\hat{V}_{n,m}^{(p,L)}$, multiplying $(p-\Delta)\hat{V}_{n,m}^{(p,L)} = 0$ by $V_k^{(p)}$, subtracting and integrating over $\Omega_L$, one gets

$$0 = \int_{\partial B_L} \left( \hat{V}_{n,m}^{(p,L)} \underbrace{\partial_n V_k^{(p)}}_{-\mathcal{M}_p^L V_k^{(p)}} - V_k^{(p)} \underbrace{\partial_n \hat{V}_{n,m}^{(p,L)}}_{=\hat{\mu}_n^{(p,L)} \hat{V}_{n,m}^{(p,L)}} \right) + \int_{\partial \Omega} \left( \hat{V}_{n,m}^{(p,L)} \underbrace{\partial_n V_k^{(p)}}_{=\mu_k^{(p)} V_k^{(p)}} - V_k^{(p)} \partial_n \hat{V}_{n,m}^{(p,L)} \right)$$

(note that the integer index $n$ should not be confusing with the normal derivative $\partial_n$). Using $\hat{V}_{n,m}^{(p,L)}|_{\partial B_L} = \psi_{n,m}$ and the self-adjointness of $\mathcal{M}_p^L$, we find the fourth identity

$$(\mu_n^{(p,L)} + \hat{\mu}_n^{(p,L)}) \int_{\partial B_L} \psi_{n,m} V_k^{(p)} = \int_{\partial \Omega} v_k^{(p)} \left( \mu_k^{(p)} \hat{V}_{n,m}^{(p,L)} - \partial_n \hat{V}_{n,m}^{(p,L)} \right). \quad (D3)$$

Note also that using the inequality$^{94}$

$$\frac{I_v(z)}{I_{v+1}(z)} < 1 + \frac{2(v+1)}{z} \quad (z > 0, v > -1) \quad (D4)$$

and setting $v = n - 2 + d/2$, we conclude that

$$\hat{\mu}_n^{(p,L)} - \hat{\mu}_n^{(0,L)} = \sqrt{p} \frac{I_{n-2+d/2}(\sqrt{p}L)}{I_{n-1+d/2}(\sqrt{p}L)} - \frac{2n+d-2}{L} \le \sqrt{p} \quad (p \ge 0) \quad (D5)$$

for any $n = 0, 1, 2, \dots$ and $d \ge 2$.

## Appendix E: Numerical method for solving the exterior Steklov problem

Various numerical methods have been developed for solving the Steklov problem in bounded domains$^{68-78}$. For instance, one can use a weak formulation of the boundary value problem (2) and basis functions on a mesh of a given bounded domain to get a matrix representation of the Dirichlet-to-Neumann operator. A numerical diagonalization of such a truncated matrix allows one to approximate a number of eigenvalues $\mu_k^{(p)}$ and to construct the associated eigenfunctions $v_k^{(p)}$ on $\partial \Omega$ and the Steklov eigenfunctions $V_k^{(p)}$ on $\Omega$. In$^{49}$, this finite-element method (FEM) was also adapted for solving mixed Steklov problems in bounded domains when the Steklov boundary condition $\partial_n V_k^{(p)} = \mu_k^{(p)} V_k^{(p)}$ on a subset $\Gamma$ of the boundary $\partial \Omega$ is completed by Dirichlet, Neumann or Robin boundary condition on the remaining part $\partial \Omega \setminus \Gamma$. In this Appendix, we discuss its extension to exterior domains.

### 1. Transparent boundary condition

We look for the solution of the exterior boundary value problem:

$$(p - \Delta)u = 0 \quad \text{in } \Omega, \quad u|_{\partial \Omega} = f, \quad u|_{\partial x \to \infty} = 0, \quad (E1)$$

for a given function $f$ on $\partial \Omega$ and $p \ge 0$. The main difficulty for implementing a FEM is that the domain $\Omega$ is unbounded.

Let $B_L$ be the ball of radius $L$, centered at the origin, that includes $\partial \Omega$. We introduce the bounded domain $\Omega_L = \Omega \cap B_L$, which lies between the inner boundary $\partial \Omega$ (with Steklov condition) and the outer boundary $\partial B_L$. As $u$ satisfies the modified Helmholtz equation in $\Omega$, it is an analytic function away from the boundary $\partial \Omega$; in particular, $u$ and its derivative are continuous on $\partial B_L$:

$$u|_{\partial B_L^-} = u|_{\partial B_L^+}, \quad (\partial_n u)|_{\partial B_L^-} = -(\partial_n u)|_{\partial B_L^+}, \quad (E2)$$

where $\partial B_L^-$ and $\partial B_L^+$ denote the approach to the sphere $\partial B_L$ from the inner and the outer sides of $\partial B_L$, respectively. By introducing an auxiliary Dirichlet-to-Neumann operator $\mathcal{M}_p^L$ in the exterior of the ball $B_L$, one gets

$$(\partial_n u)|_{\partial B_L^+} = \mathcal{M}_p^L u|_{\partial B_L^+} \Rightarrow (\partial_n u)|_{\partial B_L^-} = -\mathcal{M}_p^L u|_{\partial B_L^-}. \quad (E3)$$
In other words, one can focus on the solution of the following boundary value problem in the bounded domain $\Omega_L$:

$$ \left\{ \begin{array}{ll} (p-\Delta)u(x) = 0 & (x \in \Omega_L), \\ u(x) = f(x) & (x \in \partial\Omega), \\ \partial_n u(x) = -\mathcal{M}_p^L u(x) & (x \in \partial B_L), \end{array} \right. \qquad (E4) $$

where we wrote $\partial B_L$ instead of $\partial B_L^-$ for brevity. In this way, the original problem in the unbounded domain $\Omega$ is reduced to an equivalent problem in the bounded domain $\Omega_L$ with a Robin-like boundary condition on $\partial B_L$, known as transparent boundary condition (TBC) (see$^{79-84}$ and references therein). In particular, one can now apply the FEM to compute the eigenvalues and eigenfunctions of the Dirichlet-to-Neumann operator $\mathcal{M}_p$ in $\Omega_L$, which are identical to those in the original domain $\Omega$. A similar approach to electromagnetic and acoustic waves employed so-called perfectly matched layer condition$^{85-89}$ (see also$^{90}$).

The choice of a ball $B_L$ as an enclosing domain has the advantage that the spectral properties of the auxiliary Dirichlet-to-Neumann operator $\mathcal{M}_p^L$ are known explicitly (see Sec. III). In fact, the action of $\mathcal{M}_p^L$ onto a function $f$ on $\partial B_L$ can be written as

$$ \mathcal{M}_p^L f = \sum_{n,m} \mu_n^{(p,L)} \psi_{n,m} \int_{\partial B_L} \psi_{n,m}^* f, \qquad (E5) $$

where $\mu_n^{(p,L)}$ and $\psi_{n,m}$ are the eigenvalues and eigenfunctions of $\mathcal{M}_p^L$.

A practical implementation of the transparent boundary condition on $\partial B_L$ is straightforward. The infinite sum in Eq. (E5) has to be truncated by considering only the eigenfunctions $\psi_{n,m}$ whose eigenvalues $\mu_n^{(p,L)}$ do not exceed a chosen threshold $\mu$. We enumerate a finite number of such eigenfunctions by a single index $j$ ranging from 1 to $J$. For instance, in two dimensions, one can fix the truncation order $n_{max}$ and then take the first $J = 2n_{max} + 1$ eigenmodes, given that $\psi_{n,m} = e^{im\theta}/\sqrt{2\pi L}$, with $m = \pm 1$ and $n = 0, 1, \dots, n_{max}$. In three dimensions, one would in general have $J = (n_{max} + 1)^2$ spherical harmonics $\psi_{n,m}$ with $n \le n_{max}$. However, it is possible to make $J$ proportional to $n_{max}$ for axisymmetric domains, as described in the next subsection. For the set $\{x_1, \dots, x_N\}$ of mesh points on $\partial B_L$, we introduce the matrix $\mathbf{V}$ of size $N \times J$ with the elements $\mathbf{V}_{i,j} = \psi_{n_{j,m}}(x_i)$, and the diagonal matrix $\mathbf{m}$ of size $J \times J$ with the elements $\mathbf{m}_{i,j} = \delta_{i,j} \mu_{n_{j}}^{(p,L)}$. If $\mathbf{f}_i = f(x_i)$ are the elements of a vector $\mathbf{f}$ representing the values of a given function $f$ at points $\{x_i\}$, one gets

$$ (\partial_n u)(x_i) = -(\mathcal{M}_p^L f)(x_i) \approx -\sum_{j=1}^{J} \mu_{n_j}^{(p,L)} \mathbf{V}_{i,j} \sum_{k=1}^{N} ds_k \mathbf{V}_{k,j}^* \mathbf{f}_k = -[\mathbf{V} \mathbf{m} \mathbf{V}^\dagger \mathbf{S} \mathbf{f}]_i, \qquad (E6) $$

where $[\mathbf{V}^\dagger]_{i,j} = \mathbf{V}_{j,i}^*$, and $\mathbf{S}_{j,j'} = \delta_{j,j'} ds_j$ are the elements of the diagonal matrix $\mathbf{S}$ of size $N \times N$ formed by the surface elements $ds_j$ (e.g., in two dimensions, $ds_j$ is the distance between two neighboring points $x_j$ and $x_{j+1}$ on $\partial B_L$). In other words, the transparent boundary condition turns out to be equivalent to a generalized Robin boundary condition when the normal derivative of $u$ is proportional to an integral kernel acting on $u$ on $\partial B_L$. Such a boundary condition was implemented in$^{49}$ through a matrix $\mathbf{Q}$. Setting $\mathbf{Q} = -\mathbf{V} \mathbf{m} \mathbf{V}^\dagger \mathbf{S}$, we can thus incorporate the transparent boundary condition on $\partial B_L$.

A potential limitation of this implementation is that the eigenvalues $\mu_n^{(p,L)}$ grow with $n$ so that the approximation (E6) might be sensitive to truncation. To resolve this issue, one can implement an equivalent condition:

$$ -(\mathcal{M}_p^L)^{-1}(\partial_n u)|_{\partial B_L} = u|_{\partial B_L}, \qquad (E7) $$

that is applicable when $\mu_0^{(p,L)} > 0$ and admits a similar matrix representation. However, for all examples considered in the paper, the representation (E6) did not manifest any numerical instability so that we did not need to resort to the alternative form (E7).

Once the Steklov eigenfunctions $V_k^{(p)}$ are constructed inside the computational domain $\Omega_L$ via the above FEM, their extension outside $\Omega_L$ can be easily performed. In fact, the restriction of $V_k^{(p)}$ onto $\partial B_L$ uniquely determines its extension to $E_L = \mathbb{R}^d \setminus \overline{B_L}$ as

$$ V_k^{(p)}(\mathbf{y}) = \sum_{n,m} V_{n,m}^{(p,L)}(\mathbf{y}) \int_{\partial B_L} \psi_{n,m}^* V_k^{(p)}(\mathbf{y}) \quad (|\mathbf{y}| \ge L). \qquad (E8) $$

As previously, the infinite sum has to be truncated to a finite number of terms, whereas the integral over the sphere $\partial B_L$ is approximated by a Riemann sum on mesh points lying on $\partial B_L$. As a consequence, the evaluation of $V_k^{(p)}$ on a set of points $\mathbf{y}_1, \dots, \mathbf{y}_M$ in $E_L$ can be expressed in a matrix form, allowing for a rapid computation.
<figure><img src="image_10.png" /><figcaption>FIG. 9. Axisymmetric domains in three dimensions. (a) the exterior $\Omega = \mathbb{R}^3 \setminus \Omega_0$ of a prolate spheroid $\Omega_0$ (in gray); (b) the exterior domain $\Omega$ is split into a bounded computational domain $\Omega_L = \Omega \cap B_L$ (in light blue) and the exterior $E_L$ of a ball $B_L$; the Steklov condition is imposed on the boundary $\partial\Omega$, while the transparent boundary condition is imposed on the sphere $\partial B_L$; (c) projection of the domain $\Omega$ onto the $xz$ plane; (d) half of this projection that corresponds to $\phi = 0$ in cylindrical coordinates $(r, z, \phi)$; one can thus associate $r = x$, employ Neumann boundary condition (N) on two vertical segments, and solve the mixed Steklov problem in the bounded computational domain $\Omega_L$.</figcaption></figure>

## 2. FEM for three-dimensional axisymmetric domains

The FEM that was originally developed in<sup>49</sup> for planar domains, is directly applicable in higher dimensions. A standard practical limitation of the FEM is that a representation of the modified Helmholtz equation on a higher-dimensional mesh would involve much larger matrices and thus require substantially higher computational resources. For this reason, we focus here on three-dimensional axisymmetric domains, for which the computation can be reduced to planar domains and is thus much simpler and faster.

In practice, we assume the rotational symmetry around the z-axis (see Fig. 9) and employ the cylindrical coordinates $(r, z, \phi)$, in which the Laplace operator reads

$$ \Delta = \frac{1}{r} \partial_r r \partial_r + \partial_z^2 + \frac{1}{r^2} \partial_\phi^2, \qquad (E9) $$

where $r = \sqrt{x^2 + y^2}$. If a given function $f$ on the boundary $\partial\Omega$ of an axisymmetric domain does not depend on $\phi$, the last term in Eq. (E9) can be dropped, and the modified Helmholtz equation can be written in a matrix form as

$$ -\nabla \cdot (c\nabla u) + pru = 0, \qquad (E10) $$

where $c$ is a $2 \times 2$ diagonal matrix whose elements are given by $r$. As a consequence, the original three-dimensional problem in $\Omega$ is reduced to solving Eq. (E10) in a planar domain that is the projection of $\Omega$ onto the $xz$ plane (i.e., by setting $\phi = 0$ and $\phi = \pi$). Moreover, we further reduce this planar computational domain by considering only its half that corresponds to $\phi = 0$. To ensure the axial symmetry, one imposes an additional Neumann boundary condition at the vertical segment at $r = x = 0$. In this way, one can directly apply the FEM, as described in<sup>49</sup>. This computation can be easily extended to axisymmetric domains in higher dimensions.

## 3. Computation of non-axisymmetric eigenfunctions

For three-dimensional domains, we only presented axisymmetric Steklov eigenfunctions which do not depend on the angular coordinate $\phi$. To get a numerical solution in the general case, one needs to solve the modified Helmholtz equation in cylindrical coordinates $(r, z, \phi)$:

$$ (\Delta - p)u = \left( \frac{1}{r} \partial_r r \partial_r + \partial_z^2 + \frac{1}{r^2} \partial_\phi^2 - p \right) u = 0. \qquad (E11) $$

The axial symmetry of the domain $\Omega$ (which is still axisymmetric) allows one to search a solution in the form $u(r, z, \phi) = e^{im\phi} U_m(r, z)$, with an integer $m$, where $U_m(r, z)$ satisfies

$$ \left( \partial_r r \partial_r + r \partial_z^2 - \frac{m^2}{r} - pr \right) U_m = 0. \qquad (E12) $$
<table><caption>TABLE I. List of the first 11 eigenvalues $\mu_k^{(0)}$ (columns 2-4) and root mean squared errors (columns 5-6) for the associated eigenfunctions $v_k^{(0)}$ for the exterior of the unit disk. Similar information for $p = 1$ is provided in columns 7-11. We used $L = 2$ and $n_{max} = 30$.</caption><thead><tr><th rowspan="3">Index k</th><th colspan="6">p = 0</th><th colspan="6">p = 1</th></tr><tr><th colspan="3">Eigenvalues</th><th colspan="3">RMSE</th><th colspan="3">Eigenvalues</th><th colspan="3">RMSE</th></tr><tr><th>Exact</th><th>h<sub>max</sub> = 0.025</th><th>h<sub>max</sub> = 0.01</th><th>h<sub>max</sub> = 0.025</th><th>h<sub>max</sub> = 0.01</th><th>Exact</th><th>h<sub>max</sub> = 0.025</th><th>h<sub>max</sub> = 0.01</th><th>h<sub>max</sub> = 0.025</th><th>h<sub>max</sub> = 0.01</th></tr></thead><tbody><tr><td>0</td><td>0</td><td>1.25 × 10<sup>-9</sup></td><td>2.78 × 10<sup>-10</sup></td><td>0.0014</td><td>0.0001</td><td>1.4296</td><td>1.4302</td><td>1.4293</td><td>0.007</td><td>0.0004</td></tr><tr><td>1</td><td>1</td><td>1.0007</td><td>0.9996</td><td>0.0022</td><td>0.0013</td><td>1.6995</td><td>1.7002</td><td>1.6992</td><td>0.0059</td><td>0.0032</td></tr><tr><td>2</td><td>1</td><td>1.0007</td><td>0.9996</td><td>0.0022</td><td>0.0013</td><td>1.6995</td><td>1.7003</td><td>1.6993</td><td>0.0059</td><td>0.0032</td></tr><tr><td>3</td><td>2</td><td>2.0017</td><td>2.0000</td><td>0.0052</td><td>0.0031</td><td>2.3704</td><td>2.3723</td><td>2.3705</td><td>0.0204</td><td>0.0032</td></tr><tr><td>4</td><td>2</td><td>2.0017</td><td>2.0000</td><td>0.0052</td><td>0.0031</td><td>2.3704</td><td>2.3723</td><td>2.3705</td><td>0.0204</td><td>0.0031</td></tr><tr><td>5</td><td>3</td><td>3.0042</td><td>3.0005</td><td>0.0088</td><td>0.0060</td><td>3.2288</td><td>3.2333</td><td>3.2294</td><td>0.0235</td><td>0.0081</td></tr><tr><td>6</td><td>3</td><td>3.0042</td><td>3.0005</td><td>0.0088</td><td>0.0060</td><td>3.2288</td><td>3.2333</td><td>3.2294</td><td>0.0235</td><td>0.0081</td></tr><tr><td>7</td><td>4</td><td>4.0090</td><td>4.0013</td><td>0.0197</td><td>0.0079</td><td>4.1605</td><td>4.1699</td><td>4.1619</td><td>0.0239</td><td>0.0050</td></tr><tr><td>8</td><td>4</td><td>4.0091</td><td>4.0013</td><td>0.0197</td><td>0.0079</td><td>4.1605</td><td>4.1700</td><td>4.1619</td><td>0.0219</td><td>0.0050</td></tr><tr><td>9</td><td>5</td><td>5.0167</td><td>5.0026</td><td>0.0197</td><td>0.0111</td><td>5.1225</td><td>5.1397</td><td>5.1252</td><td>0.0248</td><td>0.0103</td></tr><tr><td>10</td><td>5</td><td>5.0167</td><td>5.0026</td><td>0.0284</td><td>0.0111</td><td>5.1225</td><td>5.1397</td><td>5.1252</td><td>0.0248</td><td>0.0103</td></tr></tbody></table>

Following the standard FEM procedure of representing the solution on the basis of functions $\varphi_j$ (see$^{49}$ for details), one gets in a matrix form

$$ U[\mathbf{K} + p\mathbf{M} + m^2\hat{\mathbf{M}}] = \mathbf{F}, \qquad (E13) $$

where the vectors $U$ and $F$ represent the solution $U_m$ and the boundary condition on basis functions, $K$ and $M$ are the usual stiff and mass matrices, whereas the matrix $\hat{\mathbf{M}}$ has the matrix elements

$$ \hat{\mathbf{M}}_{j,k} = \int_{\Omega} \frac{1}{r} \varphi_j \varphi_k. \qquad (E14) $$

In other words, one needs to replace the matrix $p\mathbf{M}$ by the matrix $p\mathbf{M} + m^2\hat{\mathbf{M}}$ in the original construction of the matrix representation of the Dirichlet-to-Neumann operator $M_p^{(49)}$. In the special case $m = 0$, the diagonalization of the matrix representing $M_p$ yields the axisymmetric eigenfunctions and the associated eigenvalues of $M_p$, as considered in the main text. In the same vein, setting $m = 1, 2, \dots$ allows one to construct other matrix representations, whereas their diagonalization will result in non-axisymmetric eigenfunctions and associated eigenvalues of $M_p$.

## 4. Validation

We validate the proposed FEM with the transparent boundary condition by solving the Steklov problem for the exterior of a disk, for which the eigenvalues are known explicitly (see Sec. III). To limit of the effect of the disk symmetry, the center of the disk is shifted from the origin. In practice, we consider a disk of radius $R = 1$ centered at $(0,0.25)$ and surrounded by a circle $\partial B_L$ of radius $L = 2$ centered at the origin. The exact eigenvalues for this exterior problem are given by $\mu_n^{(p)} = -\sqrt{p} \frac{K_n'(R\sqrt{p})}{K_n(R\sqrt{p})}$, which reduce to $\mu_n^{(0)} = n/R$ at $p = 0$. The rotational symmetry of the domain implies that the eigenfunctions are the Fourier harmonics $e^{im\theta}/\sqrt{2\pi R}$ (with $m = \pm 1$), which do not depend on $p$. Table I compares the first 11 eigenvalues to those obtained by the FEM on two meshes with different maximal mesh sizes $h_{max}$, for $p = 0$ and $p = 1$. The root mean squared errors of the associated eigenfunctions are also shown. In turn, Fig. 10 illustrates how fast the numerical eigenvalues converge to the theoretical value as the truncation order $n_{max}$ increases.

In three dimensions, we consider a sphere of radius $R = 1$ centered at $(0,0,0.25)$, surrounded by a sphere of radius $L = 2$ centered at the origin. The comparison between exact and numerically obtained eigenvalues and eigenfunctions was very similar and thus not shown.

<footnote n="1">¹J. W. S. Rayleigh, *The Theory of Sound*, 2nd ed., Vols. 1 and 2 (Dover, New York, 1945).</footnote>

<footnote n="2">²J. D. Jackson, *Classical Electrodynamics*, 3rd ed. (John Wiley & Sons, New York, 1999).</footnote>

<footnote n="3">³E. B. Davies, *Spectral Theory and Differential Operators* (Cambridge University Press, Cambridge, UK, 1995).</footnote>

<footnote n="4">⁴R. F. Bass, *Diffusions and Elliptic Operators* (Springer, New York, 1998).</footnote>

<footnote n="5">⁵C. W. Gardiner, *Handbook of stochastic methods for physics, chemistry and the natural sciences* (Springer: Berlin, 1985).</footnote>
<figure><img src="image_11.png" /><figcaption>FIG. 10. Ratio between numerical and analytical eigenvalues $\mu_1^{(0)}$ and $\mu_3^{(0)}$ for the exterior of the unit disk as a function of $n_{max}$, with $h_{max} = 0.025$.</figcaption></figure>

<footnote n="6">⁶I. N. Sneddon, *Mixed Boundary Value Problems in Potential Theory* (Wiley, New York, 1966).</footnote>

<footnote n="7">⁷Z. Schuss, *Brownian Dynamics at Boundaries and Interfaces in Physics, Chemistry and Biology* (Springer, New York, 2013).</footnote>

<footnote n="8">⁸S. Redner, *A Guide to First Passage Processes* (Cambridge: Cambridge University press, 2001).</footnote>

<footnote n="9">⁹R. Metzler, G. Oshanin, and S. Redner (Eds.) *First-Passage Phenomena and Their Applications* (Singapore: World Scientific, 2014).</footnote>

<footnote n="10">¹⁰K. Lindenberg, R. Metzler, and G. Oshanin (Eds.) *Chemical Kinetics: Beyond the Textbook* (New Jersey: World Scientific, 2019).</footnote>

<footnote n="11">¹¹D. S. Grebenkov, R. Metzler, and G. Oshanin (Eds), *Target Search Problems* (Springer: Cham, Switzerland, 2024).</footnote>

<footnote n="12">¹²D. S. Grebenkov and B.-T. Nguyen, *Geometrical structure of Laplacian eigenfunctions*, SIAM Rev. **55**, 601-667 (2013).</footnote>

<footnote n="13">¹³W. Steklov (V. A. Steklov), Sur les problèmes fondamentaux de la physique mathématique, Ann. Sci. Ec. Norm. Supér. **19**, 455-490 (1902).</footnote>

<footnote n="14">¹⁴N. Kuznetsov, T. Kulczycki, M. Kwaśnicki, A. Nazarov, S. Poborchi, I. Polterovich, and B. Siudeja, The legacy of Vladimir Andreevich Steklov, Not. Am. Math. Soc. **61**, 9-22 (2014).</footnote>

<footnote n="15">¹⁵M. Levitin, D. Mangoubi, and I. Polterovich, *Topics in Spectral Geometry* (Graduate Studies in Mathematics, vol. 237; American Mathematical Society, 2023).</footnote>

<footnote n="16">¹⁶A. Girouard and I. Polterovich, *Spectral geometry of the Steklov problem*, J. Spectr. Th. **7**, 321-359 (2017).</footnote>

<footnote n="17">¹⁷B. Colbois, A. Girouard, C. Gordon, and D. Sher, Some recent developments on the Steklov eigenvalue problem, Rev. Mat. Complut. **37**, 1-161 (2024).</footnote>

<footnote n="18">¹⁸M. Cheney, D. Isaacson, and J. C. Newell, *Electrical Impedance Tomography*, SIAM Rev. **41**, 85-101 (1999).</footnote>

<footnote n="19">¹⁹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, Rio de Janeiro, (1980), 65-73; Reprinted in Comput. Appl. Math. **25**, 2-3, 133-138 (2006).</footnote>

<footnote n="20">²⁰L. Borcea, *Electrical impedance tomography*, Inv. Prob. **18**, R99 (2002).</footnote>

<footnote n="21">²¹J. Sylvester and G. Uhlmann, *A Global Uniqueness Theorem for an Inverse Boundary Value Problem*, Ann. Math. **125**, 153-169 (1987).</footnote>

<footnote n="22">²²E. Curtis and J. Morrow, *The Dirichlet to Neumann map for a resistor network*, SIAM J. Appl. Math. **51**, 1011-1029 (1991).</footnote>

<footnote n="23">²³D. W. Fox and J. R. Kuttler, *Sloshing frequencies*, Z. Angew. Math. Phys. **34**, 668-696 (1983).</footnote>

<footnote n="24">²⁴V. Kozlov and N. Kuznetsov, *The ice-fishing problem: the fundamental sloshing frequency versus geometry of holes*, Math. Methods Appl. Sci. **27**, 289-312 (2004).</footnote>

<footnote n="25">²⁵M. Levitin, L. Parnovski, I. Polterovich, and D. A. Sher, *Sloshing, Steklov and corners: Asymptotics of Steklov eigenvalues for curvilinear polygons*, Proc. London Math. Soc. **125**, 359-487 (2022).</footnote>

<footnote n="26">²⁶D. S. Grebenkov, *Spectral theory of imperfect diffusion-controlled reactions on heterogeneous catalytic surfaces*, J. Chem. Phys. **151**, 104108 (2019).</footnote>

<footnote n="27">²⁷D. S. Grebenkov, *Probability distribution of the boundary local time of reflected Brownian motion in Euclidean domains*, Phys. Rev. E **100**, 062110 (2019).</footnote>

<footnote n="28">²⁸D. S. Grebenkov, *Paradigm shift in diffusion-mediated surface phenomena*, Phys. Rev. Lett. **125**, 078102 (2020).</footnote>

<footnote n="29">²⁹D. S. Grebenkov, *Diffusion-Controlled Reactions: An Overview, Molecules **28**, 7570 (2023).</footnote>

<footnote n="30">³⁰G. Auchmuty, Steklov eigenproblems and the representation of solutions of elliptic boundary value problems, Numer. Funct. Anal. Optim. **25**, 321-348 (2005).</footnote>

<footnote n="31">³¹G. Auchmuty and Q. Han, *Spectral representations of solutions of linear elliptic equations on exterior regions*, J. Math. Anal. Appl. **398**, 1-10 (2013).</footnote>

<footnote n="32">³²G. Auchmuty and M. Cho, *Boundary integrals and approximations of harmonic functions*, Numer. Funct. Anal. Optim. **36**, 687-703 (2015).</footnote>

<footnote n="33">³³G. Auchmuty and Q. Han, *Representations of Solutions of Laplacian Boundary Value Problems on Exterior Regions*, Appl. Math. Optim. **69**, 21-45 (2014).</footnote>

<footnote n="34">³⁴G. Auchmuty, Steklov Representations of Green's Functions for Laplacian Boundary Value Problems, Appl. Math. Optim. **77** (2018).</footnote>

<footnote n="35">³⁵R. Samson and J. M. Deutch, *Exact solution for the diffusion controlled rate into a pair of reacting sinks*, J. Chem. Phys. **67**, 847 (1977).</footnote>

<footnote n="36">³⁶S. Torquato, *Concentration dependence of diffusion-controlled reactions among static reactive sinks*, J. Chem. Phys. **85**, 7178 (1986).</footnote>

<footnote n="37">³⁷S. B. Lee, I. C. Kim, C. A. Miller, and S. Torquato, *Random-walk simulation of diffusion-controlled processes among static traps*, Phys. Rev. B **39**, 11833 (1989).</footnote>

<footnote n="38">³⁸S. Traytak, *The diffusive interaction in diffusion-limited reactions: the time-dependent case*, Chem. Phys. **193**, 351-366 (1995).</footnote>

<footnote n="39">³⁹M. Galanti, D. Fanelli, S. Traytak, and F. Piazza, *Theory of diffusion-influenced reactions in complex geometries*, Phys. Chem. Chem. Phys. **18**, 15950-15954 (2016).</footnote>

<footnote n="40">⁴⁰D. S. Grebenkov, *Diffusion toward non-overlapping partially reactive spherical traps: fresh insights onto classic problems*, J. Chem. Phys. **152**, 244108 (2020).</footnote>

<footnote n="41">⁴¹D. S. Grebenkov, *Statistics of boundary encounters by a particle diffusing outside a compact planar domain*, J. Phys. A.: Math. Theor. **54**, 015003 (2021).</footnote>

<footnote n="42">⁴²C. Xiong, *Sharp bounds for the first two eigenvalues of an exterior Steklov eigenvalue problem*, ArXiv:2304.11297v1 (2023).</footnote>

<footnote n="43">⁴³W. Arendt and A. F. M. ter Elst, *The Dirichlet-to-Neumann Operator on Exterior Domains*, Potential Anal. **43**, 313-340 (2015).</footnote>

<footnote n="45">⁴⁵L. Bundrock, A. Girouard, D. S. Grebenkov, M. Levitin, and I. Polterovich, *Exterior Steklov problems (in preparation)*.</footnote>

<footnote n="45">⁴⁵T. J. Christiansen and K. Datchev, *Low energy scattering asymptotics for planar obstacles*, Pure Appl. Anal. **5**, 767-794 (2023).</footnote>

<footnote n="46">⁴⁶L. Friedlander, *Some inequalities between Dirichlet and Neumann eigenvalues*, Arch. Rational Mech. Anal. **116**, 153-160 (1991).</footnote>

<footnote n="47">⁴⁷M. Abramowitz and I. A. Stegun, *Handbook of Mathematical Functions* (Dover Publisher, New York, 1965).</footnote>

<footnote n="48">⁴⁸D. S. Grebenkov, *Spectral properties of the Dirichlet-to-Neumann operator for spheroids*, Phys. Rev. E **109**, 055306 (2024).</footnote>
<footnote n="49">⁴⁹A. Chaigneau and D. S. Grebenkov, A numerical study of the Dirichlet-to-Neumann operator in planar domains, J. Phys. A: Math. Theor. **57**, 445201 (2024).</footnote>

<footnote n="50">⁵⁰N. Landkof, *Foundations of Modern Potential Theory* (Springer Verlag, Berlin, 1972).</footnote>

<footnote n="51">⁵¹T. Ransford and J. Rostand, Computation of capacity, Math. Comp. **76**, 1499-1520 (2007).</footnote>

<footnote n="52">⁵²T. Ransford, Computation of Logarithmic Capacity. Comput. Methods Func. Theory **10**, 555-578 (2010).</footnote>

<footnote n="53">⁵³W. Arendt and R. Mazzeo, Spectral properties of the Dirichlet-to-Neumann operator on Lipschitz domains, Ulmer Seminare **12**, 23-37 (2007).</footnote>

<footnote n="54">⁵⁴W. Arendt and R. Mazzeo, Friedlander's eigenvalue inequalities and the Dirichlet-to-Neumann semigroup, Comm. Pure Appl. Anal. **11**, 2201-2212 (2012).</footnote>

<footnote n="55">⁵⁵A. V. Skorokhod, Stochastic equations for diffusion processes in a bounded region, Theory Probab. Appl. **6**, 264 (1961).</footnote>

<footnote n="56">⁵⁶A. V. Skorokhod, Stochastic equations for diffusion processes in a bounded region. II, Theory Probab. Appl. **7**, 3 (1962).</footnote>

<footnote n="57">⁵⁷K. Itô and H. P. McKean, *Diffusion Processes and their Sample Paths*, Reprint of the 1974 edition (Springer Science & Business Media, 1996).</footnote>

<footnote n="58">⁵⁸M. I. Freidlin, *Functional Integration and Partial Differential Equations* (Princeton University Press, 1985).</footnote>

<footnote n="59">⁵⁹Y. Saisho, Stochastic differential equations for multidimensional domain with reflecting boundary, Probab. Theory Relat. Fields **74**, 455 (1987).</footnote>

<footnote n="60">⁶⁰D. S. Grebenkov, Joint distribution of multiple boundary local times and related first-passage time problems, J. Stat. Mech. 103205 (2020).</footnote>

<footnote n="61">⁶¹D. S. Grebenkov, Surface Hopping Propagator: An Alternative Approach to Diffusion-Influenced Reactions, Phys. Rev. E **102**, 032125 (2020).</footnote>

<footnote n="62">⁶²D. S. Grebenkov, An encounter-based approach for restricted diffusion with a gradient drift, J. Phys. A: Math. Theor. **55**, 045203 (2022).</footnote>

<footnote n="63">⁶³D. S. Grebenkov, Depletion of Resources by a Population of Diffusing Species, Phys. Rev. E **105**, 054402 (2022).</footnote>

<footnote n="64">⁶⁴D. S. Grebenkov, Statistics of diffusive encounters with a small target: Three complementary approaches, J. Stat. Mech. 083205 (2022).</footnote>

<footnote n="65">⁶⁵F. C. Collins and G. E. Kimball, Diffusion-controlled reaction rates, J. Coll. Sci. **4**, 425 (1949).</footnote>

<footnote n="66">⁶⁶D. S. Grebenkov, R. Metzler, and G. Oshanin, Strong defocusing of molecular reaction times results from an interplay of geometry and reaction control, Commun. Chem. **1**, 96 (2018).</footnote>

<footnote n="67">⁶⁷M. Levitin and L. Parnovski, On the principal eigenvalue of a Robin problem with a large parameter, Math. Nachrichten **281**, 272-281 (2008).</footnote>

<footnote n="68">⁶⁸A. B. Andreev and T. D. Todorov, Isoparametric finite-element approximation of a Steklov eigenvalue problem, IMA J. Numer. Anal. **24**, 309-322 (2004).</footnote>

<footnote n="69">⁶⁹Y. Yang, Q. Li and S. Li, Nonconforming finite element approximations of the Steklov eigenvalue problem, Appl. Numer. Math. **59**, 2388-2401 (2009).</footnote>

<footnote n="70">⁷⁰H. Bi and 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).</footnote>

<footnote n="71">⁷¹Q. Li and Y. Yang, A two-grid discretization scheme for the Steklov eigenvalue problem, J. Appl. Math. Comput. **36**, 129-139 (2011).</footnote>

<footnote n="72">⁷²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).</footnote>

<footnote n="73">⁷³H. Xie, A type of multilevel method for the Steklov eigenvalue problem, IMA J. Numer. Anal. **34**, 592-608 (2013).</footnote>

<footnote n="74">⁷⁴H. Bi, H. Li, and Y. Yang, An adaptive algorithm based on the shifted inverse iteration for the Steklov eigenvalue problem, Appl. Numer. Math. **105**, 64-81 (2016).</footnote>

<footnote n="75">⁷⁵E. Akhmetgaliyev, C.-Y. Kao, and B. Osting, Computational methods for extremal Steklov problems, SIAM J. Control Optim. **55**, 1226-1240 (2017).</footnote>

<footnote n="76">⁷⁶W. Alhejaili and C.-Y. Kao, Numerical studies of the Steklov eigenvalue problem via conformal mappings, Appl. Math. Comput. **347**, 785-802 (2019).</footnote>

<footnote n="77">⁷⁷O. P. Bruno and J. Galkowski, Domains Without Dense Steklov Nodal Sets, J. Fourir Anal. Appl. **26**, 45 (2020).</footnote>

<footnote n="78">⁷⁸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, Engrg. Anal. Bound. Elem. **114**, 136-147 (2020).</footnote>

<footnote n="79">⁷⁹J. B. Keller and D. Givoli, Exact non-reflecting boundary conditions, J. Comput. Phys. **82**, 172-192 (1989).</footnote>

<footnote n="80">⁸⁰D. Givoli, Non-reflecting boundary conditions, J. Comput. Phys. **94**, 1-29 (1991).</footnote>

<footnote n="81">⁸¹T. Hagstrom, Radiation boundary conditions for the numerical simulation of waves, Acta Numerica **8**, 47-106 (1999).</footnote>

<footnote n="82">⁸²D. Givoli, High-order local non-reflecting boundary conditions: a review, Wave Motion **39**, 319-326 (2004).</footnote>

<footnote n="83">⁸³Y. Achdou, C. Sabot, and N. Tchou, Transparent boundary conditions for the Helmholtz equation in some ramified domains with a fractal boundary, J. Comput. Phys. **220**, 712-739 (2007).</footnote>

<footnote n="84">⁸⁴X. Antoine, A. Arnold, C. Besse, M. Ehrhardt, and A. Schädle, A Review of Transparent and Artificial Boundary Conditions Techniques for Linear and Nonlinear Schrödinger Equations, Commun. Comput. Phys. **4**, 729-796 (2008).</footnote>

<footnote n="85">⁸⁵J. Bérenger, A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. **114**, 185-200 (1994).</footnote>

<footnote n="86">⁸⁶F. Collino and P. Monk, The perfectly matched layer in curvilinear coordinates, SIAM J. Sci. Comput. **19**, 2061-2090 (1998).</footnote>

<footnote n="87">⁸⁷E. Becache, A.-S. Bonnet-Ben Dhia, and G. Legendre, Perfectly matched layers for the convected Helmholtz equation, SIAM J. Numer. Anal. **42**, 409-433 (2004).</footnote>

<footnote n="88">⁸⁸E. Becache, A.-S. Bonnet-Ben Dhia, and G. Legendre, Perfectly matched layers for time-harmonic acoustics in the presence of a uniform flow, SIAM J. Numer. Anal. **44**, 1191-1217 (2006).</footnote>

<footnote n="89">⁸⁹V. Kalvin, Limiting absorption principle and perfectly matched layer method for Dirichlet Laplacian in quasi-cylindrical domains, SIAM J. Math. Anal. **44**, 355-382 (2012).</footnote>

<footnote n="90">⁹⁰M. Levitin and M. Marletta, A simple method of calculating eigenvalues and resonances in domains with infinite regular ends, Proc. Roy. Soc. Edinburgh A **138**, 1043-1065 (2008).</footnote>

<footnote n="91">⁹¹T. Kato, *Perturbation Theory for Linear Operators* (2nd Ed., Springer-Verlag, Berlin, 1980).</footnote>

<footnote n="92">⁹²J. Behrndt and A. F. M. ter Elst, Dirichlet-to-Neumann maps on bounded Lipschitz domains, J. Diff. Eq. **259**, 5903-5926 (2015).</footnote>

<footnote n="93">⁹³D. Bucur, P. Freitas, and J. Kennedy, The Robin problem, In *Shape optimization and spectral theory* (De Gruyter Open, Warsaw, 2017), pp 78-119.</footnote>

<footnote n="94">⁹⁴E. K. Ifantis and P. D. Siafarikas, Inequalities involving Bessel and modified Bessel functions, J. Math. Anal. Appl. **147**, 214-227 (1990).</footnote>