# A numerical study of the generalized Steklov problem in planar domains

Adrien Chaigneau, Denis S Grebenkov

► To cite this version:

Adrien Chaigneau, Denis S Grebenkov. A numerical study of the generalized Steklov problem in planar domains. Journal of Physics A: Mathematical and Theoretical, 2024, 57 (44), pp.445201. 10.1088/1751-8121/ad7fab hal-05344210

HAL Id: hal-05344210

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

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
A numerical study of the Dirichlet-to-Neumann
operator in planar domains

Adrien Chaigneau¹,*
Denis S. Grebenkov¹,†

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

November 6, 2023

Abstract

We numerically investigate the generalized Steklov problem for the modified Helmholtz equation and focus on the relation between its spectrum and the geomet- ric structure of the domain. We address three distinct aspects: (i) the asymptotic behavior of eigenvalues for polygonal domains; (ii) the dependence of the integrals of eigenfunctions on the domain symmetries; and (iii) the localization and exponen- tial decay of Steklov eigenfunctions away from the boundary for smooth shapes and in the presence of corners. For this purpose, we implemented two complementary numerical methods to compute the eigenvalues and eigenfunctions of the associated Dirichlet-to-Neumann operator for various simply-connected planar domains. We also discuss applications of the obtained results in the theory of diffusion-controlled reactions and formulate several conjectures with relevance in spectral geometry.

# 1 Introduction

The Dirichlet-to-Neumann operator [1, 2, 3, 4] plays a prominent role in applied math-
ematics, physics, engineering and medicine. One of its most common applications is
related to medical imaging and electrical impedance tomography [5, 6], in which the
electric conductivity in the bulk has to be determined from electric measurements on
the boundary and allows for instance lung function assessment. A similar technique
was used in geophysics for imaging sub-surface structures [7, 8]. A recently developed
theoretical description of diffusion-controlled reactions relies on the eigenbasis of the
Dirichlet-to-Neumann operator to decompose the underlying propagators [9, 10, 11].
Although the spectral properties of the Dirichlet-to-Neumann operator have been in-
tensively studied over the past century, there are still many open questions and un-
solved problems that explains a rapid development of this topic during the last years
[12, 13, 14, 15, 3, 16, 17, 18, 19, 20, 21].

<footnote>*adrien.chaigneau@polytechnique.edu</footnote>

<footnote>†denis.grebenkov@polytechnique.edu (Corresponding author)</footnote>
In this study, we focus on simply-connected planar bounded domains $\Omega \in \mathbb{R}^2$ with a Lipschitz boundary $\partial\Omega$. The Dirichlet-to-Neumann operator $M_p$ associates to a function $f$ on the boundary $\partial\Omega$ another function on that boundary:

$$ M_p : H^{1/2}(\partial\Omega) \to H^{-1/2}(\partial\Omega) \\ f \mapsto (\partial_n u)|_{\partial\Omega}, \qquad (1) $$

where $\partial_n$ is the normal derivative oriented outward the domain and $u(x)$ is the solution of the boundary value problem,

$$ \left\{ \begin{array}{ll} (p - \Delta)u(x) = 0 & (x \in \Omega), \\ u(x) = f(x) & (x \in \partial\Omega), \end{array} \right. \qquad (2) $$

in the Sobolev space

$$ H^1(\Omega) = \{u \in L^2(\Omega) \mid \partial_x u \in L^2(\Omega), \partial_y u \in L^2(\Omega)\}, \qquad (3) $$

where $\Delta = \partial_x^2 + \partial_y^2$ is the Laplace operator, $p \in \mathbb{R}$ is a fixed parameter and $L^2(\Omega)$ is the space of measurable and square-integrable functions on $\Omega$. The functional space $H^{1/2}(\partial\Omega)$ is the trace of $H^1(\Omega)$:

$$ \begin{aligned} H^{1/2}(\partial\Omega) &:= \operatorname{tr}(H^1(\Omega)) \\ &:= \{v \in L^2(\partial\Omega) \mid \exists u \in H^1(\Omega) : \operatorname{tr}(u) = u|_{\partial\Omega} = v\}, \end{aligned} \qquad (4) $$

and the space $H^{-1/2}(\partial\Omega)$ is the dual of $H^{1/2}(\partial\Omega)$ (see details on functional spaces in [22, 4]). In the context of diffusion-controlled reactions, the function $f$ can be thought of as a source of molecules on the boundary $\partial\Omega$, so that $M_p f$ gives their flux density on that boundary. Throughout the paper, we focus on $p \ge 0$ so that $M_p$ is a selfadjoint operator that has a discrete spectrum [4], with the eigenvalues $\mu_k^{(p)}$ and eigenfunctions $v_k^{(p)}$ satisfying

$$ M_p v_k^{(p)} = \mu_k^{(p)} v_k^{(p)} \quad (k = 0, 1, 2, \dots). \qquad (5) $$

The eigenvalues are nonnegative and the eigenfunctions form an orthonormal basis of $L^2(\partial\Omega)$. We enumerate the eigenvalues in increasing order

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

The spectrum of $M_p$ is closely related to the spectrum of the (generalized) Steklov problem [23, 24]:

$$ \left\{ \begin{array}{ll} (p - \Delta) V_k^{(p)} = 0 & (x \in \Omega), \\ \partial_n V_k^{(p)} = \mu_k^{(p)} V_k^{(p)} & (x \in \partial\Omega), \end{array} \right. \qquad (7) $$

where $V_k^{(p)}$ are the Steklov eigenfunctions. One sees that $v_k^{(p)}$ is the restriction of $V_k^{(p)}$ on $\partial\Omega$, whereas $V_k^{(p)}$ can be obtained as the unique extension of $v_k^{(p)}$ into $\Omega$:

$$ \left\{ \begin{array}{ll} (p - \Delta) V_k^{(p)} = 0 & (x \in \Omega), \\ V_k^{(p)} = v_k^{(p)} & (x \in \partial\Omega). \end{array} \right. \qquad (8) $$
In this paper, we address three questions about the spectral properties of the Dirichlet-to-Neumann operator. First, we investigate the asymptotic behavior of the eigenvalues $\mu_k^{(p)}$ as $p \to \infty$. According to [20], one has

$$ \mu_k^{(p)} \approx \sqrt{p} \quad (p \gg 1), \qquad (9) $$

for all $k$ and all bounded domains with a smooth $C^1$ boundary. The symbol $\approx$ denotes the asymptotic behavior of $\mu_k^{(p)}$ when $p$ goes to infinity; however, it also emphasizes that the left-hand side is close to the right-hand side when $p$ is large enough. In turn, in the presence of corners one can expect

$$ \mu_k^{(p)} \approx c_k \sqrt{p} \quad (p \gg 1), \qquad (10) $$

with unknown coefficients $c_k$. In Sec. 3, we reveal how $c_k$ depend on the geometry of a polygonal domain.

Second, we look at the impact of the domain symmetry onto the coefficients

$$ A_k^{(p)} = \frac{1}{\sqrt{|\partial\Omega|}} \int_{\partial\Omega} v_k^{(p)} d\boldsymbol{x}, \qquad (11) $$

where $|\partial\Omega|$ is the Lebesgue measure of $\partial\Omega$. Note that, integrating eq. (7) over $\Omega$ and using the Green's formula, one can also represent $A_k^{(p)}$ as

$$ A_k^{(p)} = \frac{p}{\mu_k^{(p)} \sqrt{|\partial\Omega|}} \int_{\Omega} V_k^{(p)} d\boldsymbol{x}. \qquad (12) $$

These coefficients play an important role in many spectral expansions (see [10] for details). When $p = 0$, one has $\mu_0^{(0)} = 0$ and $v_0^{(0)} = 1/\sqrt{|\partial\Omega|}$ that implies $A_k^{(0)} = \delta_{k,0}$ due to orthogonality of eigenfunctions $v_k^{(0)}$ to $v_0^{(0)}$, where $\delta$ is the Kronecker symbol. Moreover, for any $p \ge 0$, one also gets $A_k^{(p)} = \delta_{k,0}$ in the case of a disk due to its rotational symmetry. In Sec. 4, we study how the coefficients $A_k^{(p)}$ depend on $p$ and $k$ for various planar shapes, and discuss implications.

Third, we analyze the behavior of the Steklov eigenfunctions $V_k^{(p)}$ away from the boundary. In [12, 15, 16, 17, 18], the localization of $V_k^{(p)}$ near the boundary $\partial\Omega$ and exponentially decaying upper bounds were shown for bounded domains with real-analytic boundary $\partial\Omega$. In Sec. 5, we inspect the exponential decay of Steklov eigenfunctions away from the boundary for smooth and polygonal domains. In particular, we highlight the role of eigenvalues $\mu_k^{(p)}$ as the decay rates.

To address these questions, we compute $\mu_k^{(p)}$, $v_k^{(p)}$ and $V_k^{(p)}$ numerically for a variety of planar domains such as ellipses, triangles, rectangles, regular polygons, Koch snowflakes and randomly generated smooth shapes (see Fig. 1).

## 2 Numerical methods

There is a very limited number of domains whose symmetries allow for separation of variables and thus lead to fully explicit formulas for the eigenvalues and eigenfunctions
<figure><img src="image_2.png" /><figcaption>Figure 1: Various planar domains considered in this study and the eigenfunction $V_4^{(1)}$ for each of them. Solid black lines indicate nodal lines.</figcaption></figure>

of $M_p$ [4, 11]. In other cases, one needs to employ numerical methods to construct the Dirichlet-to-Neumann operator. For instance, one can use finite-difference or finite ele-
ment methods to discretize the problem on a regular lattice or a mesh and to construct a matrix representing the Dirichlet-to-Neumann operator $M_p$ that needs to be diagonalized to approximate the eigenvalues $\mu_k^{(p)}$ and eigenfunctions $v_k^{(p)}$. Flexibility is one of the advantages of these techniques that can deal in the same way with more general second-order elliptic operators. In turn, the need for mesh construction and large sizes of the matrices to be diagonalized are usual drawbacks. Various improvements have been proposed to overcome these limitations: an isoparametric variant of the finite element method for solving Steklov eigenvalue problems in $\mathbb{R}^d$ for second-order, self-adjoint, elliptic differential operators [25], a two-grid discretization scheme [26, 27], a finite element multi-scale discretization with an adaptive algorithm based on the shifted inverse iteration [28], an iterative multilevel approach [29], a nonconforming finite element methods [30, 31]. Another possibility is the reformulation of the Steklov eigenvalue problem in terms of an equivalent boundary integral equation [32, 33, 34], or the method of fundamental solutions for solving the Steklov and related spectral problems for the Laplace operator [35, 36]. In the planar case, one can also employ conformal mapping to transform the original domain into a simpler domain (e.g., a disk), at the price of dealing with generalized Robin boundary condition [37].

For the purpose of our study, the basic finite element method provided a sufficient accuracy and moderate computational cost. Its practical implementation is detailed in Sec. 2.1. In addition, we discuss an alternative technique based on the restriction of Green's functions (Sec. 2.2).

## 2.1 Finite element method

The first numerical method aims at representing the Dirichlet-to-Neumann operator by a matrix in two steps. First, we construct the vector $U$ representing the solution $u$ of the modified Helmholtz problem eq. (2) with a finite element method. We discretize $\Omega$ into a triangular mesh, which has $N_i$ nodes inside the domain $\Omega$, and $N_e$ nodes on the boundary $\partial\Omega$. The interior nodes are enumerated by $i = 1, \dots, N_i$, while the boundary nodes are enumerated by $i = N_i + 1, \dots, N_i + N_e$. This discretization is equivalent to projection of the weak form of the equation onto a finite-dimensional subspace of dimension $N_p = N_i + N_e$. Let $\{\phi_i\}$, with $i = 1, \dots, N_p$, be piecewise polynomial basis functions of a subspace of $\mathcal{H}^1(\Omega)$. In our implementation, each $\phi_i$ is a “hat” function that is linear on each element and takes the value 0 at all nodes $x_j$ except for $x_i$ at which it is equal to 1. This property ensures that the functions $\{\phi_i\}_{i=1,\dots,N_i}$ vanish on the boundary and thus can serve as a basis of $\mathcal{H}_0^1(\Omega)$.

Multiplying the modified Helmholtz equation $(p-\Delta)u = 0$ by a test function $\phi_j$ with $j = 1, 2, \dots, N_i$, integrating over $\Omega$, and using the Green's formula, one gets

$$ \forall j = 1, 2, \dots, N_i, \quad \int_{\Omega} p u \phi_j + \int_{\Omega} \nabla u \nabla \phi_j = 0, \qquad (13) $$

without boundary terms because $\phi_j|_{\partial\Omega} = 0$. Next, one can approximate the solution $u$ as a linear combination of basis functions:

$$ u = \sum_{i=1}^{N_p} U_i \phi_i, \qquad (14) $$
with unknown coefficients $U_i$. We split the sum as

$$ u = \sum_{i=1}^{N_i} U_i \phi_i + \underbrace{\sum_{i=N_i+1}^{N_i+N_e} F_i \phi_i}_{=f}, \qquad (15) $$

where the first term represents $u$ inside $\Omega$ and the second one incorporates the Dirichlet boundary condition $u = f$ by setting

$$ \forall i = N_i + 1, \dots, N_i + N_e, \quad U_i = F_i, \qquad (16) $$

with the coefficients $F_i$ representing $f$. Substitution of this expansion into eq. (13) yields for all $j = 1, 2, \dots, N_i$

$$ \sum_{i=1}^{N_i} \left( p \int_{\Omega} \phi_i \phi_j + \int_{\Omega} \nabla \phi_i \nabla \phi_j \right) U_i = - \sum_{i=N_i+1}^{N_i+N_e} \left( p \int_{\Omega} \phi_i \phi_j + \int_{\Omega} \nabla \phi_i \nabla \phi_j \right) F_i. \qquad (17) $$

We denote by $K$ the stiff matrix (of size $N_p \times N_p$) given by $\int_{\Omega} \nabla \phi_i \nabla \phi_j$ and by $M$ the mass matrix (of size $N_p \times N_p$) given by $\int_{\Omega} \phi_i \phi_j$. Since the basis functions $\phi_i$ vanish on all the elements that do not contain the node $x_i$, $K_{ij}$ and $M_{ij}$ are zero except if $x_i$ and $x_j$ are the vertices of the same element and thus the matrices $K$ and $M$ are very sparse. We get the matrix formulation:

$$ (pM + K)^{ii} U^i = -(pM + K)^{ie} F, \qquad (18) $$

where $(pM + K)^{ii}$ is the $N_i \times N_i$ submatrix, $(pM + K)^{ie}$ is the $N_i \times N_e$ submatrix, $F$ is the $N_e \times 1$ vector, and $U^i$ is the $N_i \times 1$ vector that can thus be found as

$$ U^i = -[(pM + K)^{ii}]^{-1} (pM + K)^{ie} F. \qquad (19) $$

Combining these coefficients with eq. (16), we construct the whole vector $U$ of coefficients $U_i$ as

$$ U = \left( \frac{-(pM + K)^{ii}}{I^{ee}} \right)^{-1} (pM + K)^{ie} F, \qquad (20) $$

where $I^{ee}$ is the identity matrix of size $N_e \times N_e$. This concludes the first step.

The second step consists in representing the action of the normal derivative. We restart from the weak formulation of $(p - \Delta)u = 0$

$$ \forall v \in \mathcal{H}^1(\Omega), \quad \int_{\Omega} puv + \int_{\Omega} \nabla u \nabla v = \int_{\partial \Omega} hv, \qquad (21) $$

with $h = \partial_n u$. Substituting eq. (14) into eq. (21), the weak formulation reads for any test function $\phi_j$ as

$$ \sum_{i=1}^{N_p} \left( p \int_{\Omega} \phi_i \phi_j + \int_{\Omega} \nabla \phi_i \nabla \phi_j \right) U_i = \sum_{i=N_i+1}^{N_p} \left( \int_{\partial \Omega} \phi_i \phi_j \right) H_i \quad (j = 1, \dots, N_i), \qquad (22) $$
where $H_i$ are the unknown coefficients representing $h$ on the basis $\{\phi_i\}$. Denoting by $M_b$ the matrix (of size $N_e \times N_e$) given by $\int_{\partial\Omega} \phi_i \phi_j$, we get the matrix formulation

$$(pM + K)^{ep} U = M_b H, \qquad (23)$$

where $(pM + K)^{ep}$ is the $N_e \times N_p$ submatrix. As a consequence, eq. (20) for $U$ implies

$$H = M_p F, \qquad (24)$$

where

$$M_p = M_b^{-1} (pM + K)^{ep} \left( - \frac{[(pM + K)^{ii}]^{-1} (pM + K)^{ie}}{I^{ee}} \right). \qquad (25)$$

According to eq. (24) the matrix $M_p$ transforms the Dirichlet boundary condition $u = f$, with a function $f$ represented by the vector $F$, into the Neumann boundary condition $\partial_n u = h$, with $h$ represented by the vector $H$. In other words, this is a matrix representation of the Dirichlet-to-Neumann operator $M_p$ in terms of basis functions $\{\phi_i\}$. Once the matrix $M_p$ is constructed, one can apply standard numerical algorithms to diagonalize it. The obtained eigenvalues of $M_p$ approximate the eigenvalues $\mu_k^{(p)}$. As $M_p$ is a finite-size matrix (of size $N_e \times N_e$), only a finite number of eigenvalues $\mu_k^{(p)}$ can be accurately approximated. In practice, we will limit our analysis to few tens of eigenvalues (say with $k$ from 0 up to 20), for which the method is very accurate, as checked below. In turn, each eigenvector $V_k$ of the matrix $M_p$ determines the coefficients of the expansion of $v_k^{(p)}$ on the basis functions:

$$v_k^{(p)}(\mathbf{x}) = \sum_{i=1}^{N_e} \mathbf{V}_{k,i} \phi_{i+N_i}(\mathbf{x}). \qquad (26)$$

As $\phi_i$ are chosen to be the hat functions, one simply has $v_k^{(p)}(x_{j+N_i}) = V_{k,j}$ at the nodes $x_{j+N_i}$ of the boundary. As eigenvectors $V_k$ and thus the eigenfunctions $v_k^{(p)}$ are defined up to a multiplicative factor, we explicitly renormalize them to ensure the unit $L^2(\partial\Omega)$ norm, i.e. $\int_{\partial\Omega} |v_k^{(p)}|^2 = 1$. The Steklov eigenfunctions $V_k^{(p)}$ are obtained by replacing $F$ by $v_k^{(p)}$ in eq. (20).

The flexibility of the finite element method allows for various extensions. For instance, one can solve the mixed Steklov problem when the boundary $\partial\Omega$ is composed of two disjoint parts: $\partial\Omega = \partial\Omega_1 \cup \partial\Omega_2$. The Dirichlet-to-Neumann operator then acts as

$$M_p : H^{1/2}(\partial\Omega_1) \to H^{-1/2}(\partial\Omega_1) \\ f \mapsto (\partial_n u)|_{\partial\Omega_1}, \qquad (27)$$

where

$$\left\{ \begin{array}{l} (p - \Delta)u(\mathbf{x}) = 0 \qquad (\mathbf{x} \in \Omega), \\ u(\mathbf{x}) = f(\mathbf{x}) \qquad (\mathbf{x} \in \partial\Omega_1), \\ u(\mathbf{x}) = 0 \qquad (\mathbf{x} \in \partial\Omega_2), \end{array} \right. \quad \text{or} \quad \left\{ \begin{array}{l} (p - \Delta)u(\mathbf{x}) = 0 \qquad (\mathbf{x} \in \Omega), \\ u(\mathbf{x}) = f(\mathbf{x}) \qquad (\mathbf{x} \in \partial\Omega_1), \\ \partial_n u(\mathbf{x}) = 0 \qquad (\mathbf{x} \in \partial\Omega_2). \end{array} \right.$$

In other words, one imposes an additional Dirichlet or Neumann boundary condition on $\partial\Omega_2$ and defines the operator $M_p$ to act on functions on $\partial\Omega_1$. In the context of diffusion-controlled reactions, mixed boundary conditions allow one to describe various processes
in which the diffusing particle can leave the confining domain through an escape region or be destroyed on it, before reaching the target region [38]. The implementation of the Dirichlet boundary condition to our numerical method is straightforward: one just needs to truncate the second term in eq. (15) to the indices corresponding to the nodes on $\partial\Omega_1$. In this way, one incorporates the Dirichlet boundary condition $u = f$ on $\partial\Omega_1$, while letting the homogeneous Dirichlet boundary condition $u = 0$ on $\partial\Omega_2$. The inclusion of the Neumann boundary condition is as well simple: one extends the first term in eq. (15) to the indices corresponding to the nodes on $\partial\Omega_2$, i.e. the first term then represents $u$ inside $\Omega$ and $u$ on the reflecting boundary $\partial\Omega_2$, while the second term (reduced to the indices corresponding to the nodes on $\partial\Omega_1$) encodes the Dirichlet boundary condition $u = f$ on $\partial\Omega_1$. An extension of this method to three-dimensional domains is also straightforward.

## 2.2 Green's function method

The second method was inspired by the spectral decompositions of the Green's functions [10]. Let us introduce the Green's function $\tilde{G}_q(\mathbf{x}, p|\mathbf{x}_0)$ satisfying

$$ \left\{ \begin{array}{l} (p - \Delta)\tilde{G}_q(\mathbf{x}, p|\mathbf{x}_0) = \delta(\mathbf{x} - \mathbf{x}_0) \quad (\mathbf{x} \in \Omega), \\ -\partial_n \tilde{G}_q(\mathbf{x}, p|\mathbf{x}_0) = q\tilde{G}_q(\mathbf{x}, p|\mathbf{x}_0) \quad (\mathbf{x} \in \partial\Omega), \end{array} \right. \qquad (28) $$

with a constant $0 \le q \le \infty$, and $\delta(\mathbf{x} - \mathbf{x}_0)$ being the Dirac distribution. Its expansion on the Steklov eigenbasis reads [10]:

$$ \tilde{G}_q(\mathbf{x}, p|\mathbf{x}_0) = \tilde{G}_{\infty}(\mathbf{x}, p|\mathbf{x}_0) + \sum_{k=0}^{\infty} \frac{[V_k^{(p)}(\mathbf{x}_0)]^* V_k^{(p)}(\mathbf{x})}{\mu_k^{(p)} + q}, \qquad (29) $$

which is also possible to write as

$$ \tilde{G}_q(\mathbf{x}, p|\mathbf{x}_0) = \tilde{G}_0(\mathbf{x}, p|\mathbf{x}_0) - \sum_{k=0}^{\infty} \frac{[V_k^{(p)}(\mathbf{x}_0)]^* V_k^{(p)}(\mathbf{x})}{(\mu_k^{(p)}/q + 1)\mu_k^{(p)}}, \qquad (30) $$

where asterisk denotes the complex conjugate. Setting $q = 0$ and restricting the points $\mathbf{x}$ and $\mathbf{x}_0$ onto the boundary $\partial\Omega$ in eq. (29), one gets the integral kernel of the inverse of the Dirichlet-to-Neumann operator:

$$ \tilde{G}_0(\mathbf{x}, p|\mathbf{x}_0) = \sum_{k=0}^{\infty} \frac{[v_k^{(p)}(\mathbf{x}_0)]^* v_k^{(p)}(\mathbf{x})}{\mu_k^{(p)}} \quad (\mathbf{x}, \mathbf{x}_0 \in \partial\Omega). \qquad (31) $$

However, the Green's function $\tilde{G}_0(\mathbf{x}, p|\mathbf{x}_0)$ exhibits a singularity at $\mathbf{x} = \mathbf{x}_0$ that would require additional regularization and may enhance numerical errors. It is therefore convenient to remove this singularity by using eq. (30). Restricting both $\mathbf{x}$ and $\mathbf{x}_0$ to the boundary, we consider the integral kernel:

$$ \tilde{g}_q(\mathbf{x}, p|\mathbf{x}_0) = \frac{\tilde{G}_0(\mathbf{x}, p|\mathbf{x}_0) - \tilde{G}_q(\mathbf{x}, p|\mathbf{x}_0)}{q} = \sum_{k=0}^{\infty} \frac{[v_k^{(p)}(\mathbf{x}_0)]^* v_k^{(p)}(\mathbf{x})}{\mu_k^{(p)}(\mu_k^{(p)} + q)} \quad (\mathbf{x}, \mathbf{x}_0 \in \Omega). $$

One gets then

$$ \int_{\partial\Omega} \tilde{g}_q(\mathbf{x}, p|\mathbf{x}_0) v_k^{(p)}(\mathbf{x}) d\mathbf{x} = \eta_k^{(p)} v_k^{(p)}(\mathbf{x}_0), \qquad (32) $$
where the eigenvalues $\eta_k^{(p)}$ of this integral operator are related to $\mu_k^{(p)}$ as

$$ \mu_k^{(p)} = \sqrt{1/\eta_k^{(p)} + q^2/4 - q/2}. \qquad (33) $$

The discretization of the integral at boundary points $x_i$ yields a system of linear equations

$$ \sum_j \delta x_j \tilde{g}_q(x_j, p|x_i) v_k^{(p)}(x_j) = \eta_k^{(p)} v_k^{(p)}(x_i), \qquad (34) $$

where $\delta x_j$ are the areas of the surface elements around $x_j$. In other words, one needs to diagonalize the matrix $G_{ij} = \delta x_j \tilde{g}_q(x_j, p|x_i)$.

Each Green's function can be found either directly (e.g., by a finite element method), or from its (truncated) spectral expansion

$$ \tilde{G}_q(x, p|x_0) = \sum_{k=1}^{\infty} \frac{[u_k^{(q)}(x_0)]^* u_k^{(q)}(x)}{p + \lambda_k^{(q)}}, \qquad (35) $$

where $\lambda_k^{(q)}$ and $u_k^{(q)}$ are the eigenvalues and the $L^2(\Omega)$-normalized eigenfunctions of the Laplace operator:

$$ \left\{ \begin{array}{ll} -\Delta u_k^{(q)} = \lambda_k^{(q)} u_k^{(q)} & (x \in \Omega), \\ -\partial_n u_k^{(q)} = q u_k^{(q)} & (x \in \partial \Omega). \end{array} \right. \qquad (36) $$

We stress that the use of the spectral expansion eq. (35) is not optimal from the numerical point of view. In fact, the construction of Laplacian eigenfunctions and eigenvalues is a time-consuming procedure, while the truncation of the infinite series in eq. (35) to a finite number of terms can be the major source of numerical errors. We used this expansion as a straightforward way to access the Green's function for validation purposes, but more efficient numerical tools can be designed for this task. In turn, the advantage of this method is that the Laplacian eigenfunctions need to be computed only once for a given domain.

Once $\mu_k^{(p)}$ and $v_k^{(p)}$ are found, one can also compute the extension $V_k^{(p)}(x)$ by using the Green's function and the Neumann boundary condition $\partial_n V_k^{(p)} = \mu_k^{(p)} v_k^{(p)}$, so that

$$ V_k^{(p)}(x_0) = \int_{\partial \Omega} \tilde{G}_0(x, p|x_0) \mu_k^{(p)} v_k^{(p)}(x) dx. \qquad (37) $$

Note that this relation is not applicable for $p = 0$ and $k = 0$, for which $\mu_0^{(0)} = 0$ and $V_0^{(0)}(x_0) = 1/\sqrt{|\partial\Omega|}$ are already known.

## 2.3 Numerical validation

Both methods were implemented in Matlab. In particular, we relied on the Matlab PDEtool to generate triangular meshes, the matrices $K$, $M$, and $M_b$ in the first method, and the Laplacian eigenfunctions in the second method. We also used build-in functions `eig` and `eigs` for matrix diagonalizations.
We validate both numerical methods by comparing their results with the explicit formulas known for the disk of radius R [11]:

$$ \mu_0^{(p)} = \sqrt{p} \frac{I_1(R\sqrt{p})}{I_0(R\sqrt{p})}, \quad v_0^{(p)}(\theta) = \frac{1}{\sqrt{2\pi R}}, \quad V_0^{(p)}(r, \theta) = \frac{I_0(r\sqrt{p})}{I_0(R\sqrt{p})} v_0^{(p)}(\theta), \quad (38a) $$

$$ \mu_{2k}^{(p)} = \sqrt{p} \frac{I'_{2k}(R\sqrt{p})}{I_{2k}(R\sqrt{p})}, \quad v_{2k}^{(p)}(\theta) = \frac{\cos(k\theta)}{\sqrt{\pi R}}, \quad V_{2k}^{(p)}(r, \theta) = \frac{I_{2k}(r\sqrt{p})}{I_{2k}(R\sqrt{p})} v_{2k}^{(p)}(\theta), \quad (38b) $$

$$ \mu_{2k+1}^{(p)} = \sqrt{p} \frac{I'_{2k+1}(R\sqrt{p})}{I_{2k+1}(R\sqrt{p})}, \quad v_{2k+1}^{(p)}(\theta) = \frac{\sin(k\theta)}{\sqrt{\pi R}}, \quad V_{2k+1}^{(p)}(r, \theta) = \frac{I_{2k+1}(r\sqrt{p})}{I_{2k+1}(R\sqrt{p})} v_{2k+1}^{(p)}(\theta), \quad (38c) $$

where $I_k(z)$ is the modified Bessel function of the first kind, prime denotes the derivative with respect to the argument, and we use polar coordinates $(r, \theta)$. Note that all the eigenvalues are twice degenerate, except for the first one $\mu_0^{(p)}$ which is simple. For Method 1 the mesh is composed of 384138 triangles, and the maximal mesh size is 0.0042. For Method 2 the mesh is composed of 16256 triangles and the series in eq. (35) was truncated to 131 eigenfunctions of the Laplace operator. Table 1 summarizes the first 11 eigenvalues of $M_p$ for $R = 1$ and $p = 1$, while Fig. 1a shows the corresponding eigenfunction $V_4^{(1)}$. One sees that the eigenvalues in the third column, which were numerically obtained by Method 1, are in excellent agreement with the exact ones given by eq. (38). While the numerical eigenvalues computed by Method 2 are less accurate, they were computed much faster on a mesh with a smaller number of triangles. We conclude that two methods provide complementary numerical tools to access the spectral properties of the Dirichlet-to-Neumann operator. As a systematic comparison of two methods is beyond the scope of the paper, we use Method 1 in the following computations. Table 2 presents the root mean squared errors (RMSE) between the exact eigenfunctions from eq. (38) and the numerical ones. For each index $k$ we compute the RMSE from the formula: $\sqrt{N^{-1} \sum_{i=1}^{N} (f_i - f_i^*)^2}$, where $N$ is the number of points on the boundary, and $f$ and $f^*$ represent the analytical and numerical values. Expectedly, this error increases with the index $k$ but remains negligible for the considered range of indices. In Appendix 7.1, we provide additional verifications by presenting the explicitly known formulas for a rectangle and comparing them with our numerical results.

In the following, we set the maximal mesh size to 0.005.

# 3 Asymptotic behavior of eigenvalues

In this section, we investigate the asymptotic behavior of the eigenvalues $\mu_k^{(p)}$ in the limit $p \to \infty$. In particular, we focus on the dependence of $\mu_k^{(p)}$ on the angles of polygonal domains.

## 3.1 Ellipses

We start the study by considering smooth anisotropic domains such as an ellipse with semiaxes $a$ and $b$ (Fig. 1b):

$$ \Omega = \{(x, y) \in \mathbb{R}^2 : (x/a)^2 + (y/b)^2 < 1\}. \quad (39) $$
<table><caption>Table 1: List of the first 11 eigenvalues $\mu_k^{(p)}$ for the unit disk, with $p = 1$.</caption><thead><tr><th>Index k</th><th>Exact</th><th>Method 1</th><th>Method 2</th></tr></thead><tbody><tr><td>0</td><td>0.4464</td><td>0.4464</td><td>0.4464</td></tr><tr><td>1</td><td>1.2402</td><td>1.2402</td><td>1.2402</td></tr><tr><td>2</td><td>1.2402</td><td>1.2402</td><td>1.2402</td></tr><tr><td>3</td><td>2.1633</td><td>2.1633</td><td>2.1640</td></tr><tr><td>4</td><td>2.1633</td><td>2.1633</td><td>2.1640</td></tr><tr><td>5</td><td>3.1235</td><td>3.1235</td><td>3.1257</td></tr><tr><td>6</td><td>3.1235</td><td>3.1235</td><td>3.1259</td></tr><tr><td>7</td><td>4.0992</td><td>4.0993</td><td>4.1059</td></tr><tr><td>8</td><td>4.0992</td><td>4.0993</td><td>4.1061</td></tr><tr><td>9</td><td>5.0828</td><td>5.0832</td><td>5.0955</td></tr><tr><td>10</td><td>5.0828</td><td>5.0832</td><td>5.0959</td></tr></tbody></table>

<table><caption>Table 2: List of root mean squared errors for the first 11 eigenfunctions $v_k^{(p)}$ of the unit disk, with $p = 1$.</caption><thead><tr><th>Index k</th><th>RMSE 1</th><th>RMSE 2</th></tr></thead><tbody><tr><td>0</td><td>0.0000</td><td>0.0000</td></tr><tr><td>1</td><td>0.0013</td><td>0.0004</td></tr><tr><td>2</td><td>0.0013</td><td>0.0004</td></tr><tr><td>3</td><td>0.0021</td><td>0.0006</td></tr><tr><td>4</td><td>0.0021</td><td>0.0055</td></tr><tr><td>5</td><td>0.0044</td><td>0.0054</td></tr><tr><td>6</td><td>0.0044</td><td>0.0054</td></tr><tr><td>7</td><td>0.0024</td><td>0.0173</td></tr><tr><td>8</td><td>0.0024</td><td>0.0173</td></tr><tr><td>9</td><td>0.0031</td><td>0.0120</td></tr><tr><td>10</td><td>0.0031</td><td>0.0120</td></tr></tbody></table>

In this configuration, we numerically check the expected asymptotic relation eq. (9) at large $p$ as well as the opposite limit $p \to 0$ [39, 40]:

$$ \mu_0^{(p)} \approx p \frac{|\Omega|}{|\partial\Omega|} \quad (p \ll 1). \qquad (40) $$

Figure 2 illustrates these relations for two ellipses of aspect ratio 2 and 10. For the second (more elongated) ellipse (Fig. 2b), the first two eigenvalues reach the asymptotic behavior $\sqrt{p}$ slower than the others, showing an extended transient regime at intermediately large $p$. The existence of this transient regime suggests that the dependence of eigenvalues on $p$ may reveal additional geometrical features of the domains, such as its anisotropy.

<figure><img src="image_3.png" /><figcaption>Figure 2: Dependence of the first eigenvalues $\mu_k^{(p)}$ (in symbol) on $p$ for an ellipse with semiaxes 1 and 0.5 (a) and an ellipse with semiaxes 1 and 0.1 (b). Solid black line presents the asymptotic relation (9) and dotted black line indicates the relation (40).</figcaption></figure>
## 3.2 Rectangles and regular polygons

Then, we inspect the role of angles of a polygonal domain on the asymptotic behavior of the eigenvalues. For this purpose, we consider rectangular and regular polygonal domains. For a square, as $p \to \infty$, we observe the asymptotic behavior eq. (10), with $c_k \approx 0.51$ for $k = \{0, 1, 2, 3\}$, and $c_k = 1$ for other $k$. The prefactor $c_k$ was estimated from the ratio $\mu_k^{(p)}/\sqrt{p}$ at $p = 10^3$. In other words, the first four eigenvalues exhibit the asymptotic behavior eq. (10) and deviate from eq. (9), which was demonstrated for bounded domains with smooth boundaries. Note that $\mu_1^{(p)}$ and $\mu_2^{(p)}$ are identical, i.e., they are degenerate eigenvalues of multiplicity 2. In the case of a regular polygon with $N$ vertices, the first $N$ eigenvalues exhibit the asymptotic behavior eq. (10), with the same coefficients $c_0 = ... = c_{N-1} < 1$. As the number of sides of the polygon increases, the domain is getting closer to a disk so that the coefficients $c_0 = ... = c_{N-1}$ approach 1 (see Fig. 3a). We conjecture that these coefficients only depend on the angle $\alpha = \pi(1 - 2/N)$ of the regular polygon and take the value

$$c_0 = ... = c_{N-1} = \sin(\alpha/2). \qquad (41)$$

This conjecture is numerically confirmed on Fig. 3a.

<figure><img src="image_4.png" /><figcaption>Figure 3: The coefficient $c_0$ (in symbol) as a function of (a) the number $N$ of sides of the regular polygon; (b) the aspect ratio of the rectangle with one side of length 1 and the other varying from 1 to 10. Solid line presents eq. (41), while dashed horizontal line indicates the value 1 of this coefficient for domains with smooth boundaries.</figcaption></figure>

To underline the role of the angles, we plot the coefficient $c_0$ as a function of the aspect ratio of a rectangle (Fig. 3b). It appears that the anisotropy of the domain does not affect the coefficient $c_0$, whose minor variations can be attributed to weak inaccuracies of the numerical method for elongated domains.

## 3.3 Prefractals

To investigate the effect of roughness of the boundary, we consider a family of prefractal Koch snowflakes. These domains are constructed iteratively, starting from the equilateral
triangle of sidelength 2 (generation 0, Fig. 1g) and adding finer geometric features at each iteration (see Figs. 1j, 1k, 1l for generations 1, 2 and 3 respectively). Figure 4a shows the dependence of the eigenvalues on $p$ for an equilateral triangle of sidelength 2. As $p \to \infty$, we observe the asymptotic behavior eq. (10), with $c_k \approx \sin(\pi/6) = 0.5$ for $k = \{0, 1, 2\}$, and $c_k = 1$ for other $k$. This is in agreement with the conjectured expression eq. (41) of $c_k$ for regular polygons. For next three generations (see Fig. 4), it appears that the first 6, 18 and 66 eigenvalues respectively have $c_k \approx \sin(\pi/6) = 0.5$ for $k$ from 0 to $K_g - 1$, and $c_k = 1$ for other $k$, where $K_g = 3(1 + 4^0 + 4^1 + \dots + 4^{g-1}) = 2 + 4^g$, is the number of angles $\pi/3$ in the generation $g$. We conclude that the value of the coefficient $c_0$ is independent of the fractal generation $g$, which only influences the number of $c_k$ that are smaller than 1. Note also that for the considered prefactal domains, the asymptotic behavior eq. (40) still holds, which suggests that the boundary $\partial\Omega$ is not rough enough to break this asymptotic relation.

<figure><img src="image_5.png" /><figcaption>Figure 4: Dependence of some eigenvalues $\mu_k^{(p)}$ (in symbol) on $p$ for the first four generations of the Koch snowflake. Solid black line presents the asymptotic relation (9) and dotted black line indicates the relation (40).</figcaption></figure>
## 3.4 Generic triangle

We investigate the asymptotic behavior of eigenvalues for a generic triangle, which is constructed by setting the length of one side to be 2 and two angles to be $\pi/12$ and $\pi/3$ (Fig. 1d). To ensure the quality of the numerical results, we refine the mesh and set the maximal mesh size to 0.003. The last column of Table 3 presents the obtained coefficients $c_k$.

<table><caption>Table 3: First 9 coefficients $c_k$ for a triangle with angles $\pi/12$, $\pi/3$ and $7\pi/12$. The last column presents the values $c_k$ estimated from the ratio $\mu_k^{(p)}/\sqrt{p}$ at $p = 10^3$, while the previous column shows conjectured values. Columns 2-4 present the effective angles $\alpha_0$, $\alpha_1$, $\alpha_2$ (those that are equal to or exceed $\pi$ are shown in gray).</caption><thead><tr><th>k</th><th>α<sub>0</sub></th><th>α<sub>1</sub></th><th>α<sub>2</sub></th><th>c<sub>k</sub>(conjecture)</th><th>c<sub>k</sub>(numeric)</th></tr></thead><tbody><tr><td>0</td><td>π/12</td><td>π/3</td><td>7π/12</td><td>sin(1/2 π/12) = 0.1305</td><td>0.1305</td></tr><tr><td>1</td><td>3π/12</td><td>π/3</td><td>7π/12</td><td>sin(1/2 3π/12) = 0.3827</td><td>0.3833</td></tr><tr><td>2</td><td>5π/12</td><td>π/3</td><td>7π/12</td><td>sin(1/2 π/3) = 0.5</td><td>0.45999</td></tr><tr><td>3</td><td>5π/12</td><td>π</td><td>7π/12</td><td>sin(1/2 5π/12) = 0.6088</td><td>0.6104</td></tr><tr><td>4</td><td>7π/12</td><td>π</td><td>7π/12</td><td>sin(1/2 7π/12) = 0.7934</td><td>0.7935</td></tr><tr><td>5</td><td>7π/12</td><td>π</td><td>7π/4</td><td>sin(1/2 7π/12) = 0.7934</td><td>0.7957</td></tr><tr><td>6</td><td>9π/12</td><td>π</td><td>7π/4</td><td>sin(1/2 9π/12) = 0.9239</td><td>0.9260</td></tr><tr><td>7</td><td>11π/12</td><td>π</td><td>7π/4</td><td>sin(1/2 11π/12) = 0.9914</td><td>0.9925</td></tr><tr><td>8</td><td>13π/12</td><td>π</td><td>7π/4</td><td>sin(π) = 1</td><td>1.0020</td></tr></tbody></table>

To interpret the numerical results, we introduce the notion of “effective angles”. Let us denote the angles of the triangle as $\alpha_0^{(0)} = \pi/12$, $\alpha_1^{(0)} = \pi/3$ and $\alpha_2^{(0)} = 7\pi/12$. We conjecture that the coefficient $c_0$ is given by the smallest angle $\alpha_0^{(0)}$: $c_0 = \sin(\alpha_0^{(0)}/2) = \sin(\pi/24) \approx 0.1305$. The next coefficient $c_1$ turns out to be close to $c_1 \approx \sin(3\alpha_0^{(0)}/2) = \sin(\pi/8) \approx 0.3827$, i.e. one replaces the original angle $\alpha_0^{(0)}$ by $3\alpha_0^{(0)}$, by adding $2\alpha_0^{(0)}$. In
other words, one now deals with effective angles $\alpha_0^{(1)} = 3\pi/12$, $\alpha_1^{(1)} = \pi/3$, $\alpha_2^{(1)} = 7\pi/12$, and the smallest of them determines $c_1$. Then, we apply this iterative re-adjustment of the effective angles, i.e. we identify the index $i$ of the smallest effective angle and update it by adding twice larger original angle with index $i$:

$$ \alpha_i^{(k+1)} = \alpha_i^{(k)} + 2\alpha_i^{(0)}. \qquad (42) $$

For instance, we get for $k = 2$: $\alpha_0^{(2)} = 5\pi/12$, $\alpha_1^{(2)} = \pi/3$, $\alpha_2^{(2)} = 7\pi/12$, so that $c_2$ is determined by the smallest angle $\alpha_1^{(2)}$: $c_2 \approx \sin(\alpha_1^{(2)}/2) = \sin(\pi/6) = 0.5$. However, when an effective angle exceeds $\pi$, it does not contribute to the iterative procedure anymore. The effective angles and the conjectured values of $c_k$ are given in Table 3.

## 3.5 Conjecture for polygonal domains

The above iterative process can be generalized to arbitrary polygonal domains. Let $a_0 = \{\alpha_0, \dots, \alpha_{N-1}\}$ be a sequence of all angles of a polygonal domain. The coefficient $c_0$ is set to be $c_0 = \sin(\min\{\pi, a_0\}/2)$, i.e., it is determined by the smallest angle of the domain, say, $\alpha_i$. After that, the $i$-th element of the sequence $a_0$ is increased by $2\alpha_i$ to get the updated sequence $a_1$ that determines the next coefficient: $c_1 = \sin(\min\{\pi, a_1\}/2)$. At step $k$, one sets $c_k = \sin(\min\{\pi, a_k\}/2)$, i.e., the smallest angle in the already constructed sequence $a_k$. If this angle stands on the position $j$, the $j$-th element of this sequence is increased by $2\alpha_j$ to produce a new sequence $a_{k+1}$, and so on. Note that if there are several equal minima in the sequence of effective angles, one choose the one, for which the increment $2\alpha_i$ is the largest (and if the increments are identical, one chooses any of them). It is clear that, after a number of steps, all angles in the sequence will exceed $\pi$, so that all the remaining coefficients $c_k$ become equal 1. In the case of a regular polygon with $N$ sides and equal angles $\alpha$, the conjecture simply implies that the first $N$ coefficients are $c_k = \sin(\alpha/2)$, while the other $c_k = 1$, in agreement with our earlier numerical results. We check the conjecture by considering a generic polygonal domain with angles $\pi/12, \pi/12, \pi/4, \pi/4, 1.9064, 2.7224, 23\pi/12, 23\pi/12$ (Fig. 1f). Table 4 presents the effective angles, the conjectured coefficients $c_k$ and those obtained numerically.

Note that the smallest effective angle that determines the coefficient $c_k$, also indicates the region in which the related eigenfunctions $V_k^{(p)}$ is concentrated in the limit $p \to \infty$. In our examples, the value $p = 10^3$ provided very good agreement between theoretical and numerical results but one may need to increase the value of $p$ for other domains.

## 4 Coefficients of spectral expansions

In this section, we investigate how the coefficients $A_k^{(p)}$ defined by eq. (11) depend on $p$ and $k$ for various planar domains. As mentioned in Sec. 1, $A_k^{(0)} = \delta_{k,0}$ due to orthogonality of eigenfunctions $v_k^{(0)}$ to $v_0^{(0)} = 1/\sqrt{|\partial\Omega|}$. However, for $p > 0$, the eigenfunction $v_0^{(p)}$ is in general not constant so that $A_k^{(p)}$ may be non zero for $k > 0$. In the case of a disk, the rotational symmetry implies that $v_k^{(p)}$ do not depend on $p$, so that $A_k^{(p)} = \delta_{k,0}$ for any $p \ge 0$. As a consequence, spectral expansions involving $A_k^{(p)}$ as coefficients are reduced to a single term [10]. One may wonder how the domain shape can alter this behavior.
<table><caption>Table 4: First 18 coefficients $c_k$ for the polygon shown in Fig. 1f with angles $\pi/12$, $\pi/12$, $\pi/4$, $\pi/4$, 1.9064, 2.7224, $23\pi/12$, $23\pi/12$. The last column presents the values $c_k$ estimated from the ratio $\mu_k^{(p)}/\sqrt{p}$ at $p = 10^3$, while the previous column shows conjectured values. Columns 2-9 present the effective angles (those that exceed $\pi$ are shown in gray).</caption><thead><tr><th>k</th><th>α<sub>0</sub></th><th>α<sub>1</sub></th><th>α<sub>2</sub></th><th>α<sub>3</sub></th><th>α<sub>4</sub></th><th>α<sub>5</sub></th><th>α<sub>6</sub></th><th>α<sub>7</sub></th><th>c<sub>k</sub>(conjecture)</th><th>c<sub>k</sub>(numeric)</th></tr></thead><tbody><tr><td>0</td><td>π/12</td><td>π/12</td><td>π/4</td><td>π/4</td><td>1.9064</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 π/12) = 0.1305</td><td>0.1306</td></tr><tr><td>1</td><td>3π/12</td><td>π/12</td><td>π/4</td><td>π/4</td><td>1.9064</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 π/12) = 0.1305</td><td>0.1306</td></tr><tr><td>2</td><td>3π/12</td><td>3π/12</td><td>π/4</td><td>π/4</td><td>1.9064</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 π/4) = 0.3827</td><td>0.3828</td></tr><tr><td>3</td><td>3π/12</td><td>3π/12</td><td>3π/4</td><td>π/4</td><td>1.9064</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 π/4) = 0.3827</td><td>0.3828</td></tr><tr><td>4</td><td>3π/12</td><td>3π/12</td><td>3π/4</td><td>3π/4</td><td>1.9064</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 3π/12) = 0.3827</td><td>0.3838</td></tr><tr><td>5</td><td>5π/12</td><td>3π/12</td><td>3π/4</td><td>3π/4</td><td>1.9064</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 3π/12) = 0.3827</td><td>0.3838</td></tr><tr><td>6</td><td>5π/12</td><td>5π/12</td><td>3π/4</td><td>3π/4</td><td>1.9064</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 5π/12) = 0.6088</td><td>0.6113</td></tr><tr><td>7</td><td>7π/12</td><td>5π/12</td><td>3π/4</td><td>3π/4</td><td>1.9064</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 5π/12) = 0.6088</td><td>0.6113</td></tr><tr><td>8</td><td>7π/12</td><td>7π/12</td><td>3π/4</td><td>3π/4</td><td>1.9064</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 7π/12) = 0.7934</td><td>0.7967</td></tr><tr><td>9</td><td>9π/12</td><td>7π/12</td><td>3π/4</td><td>3π/4</td><td>1.9064</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 7π/12) = 0.7934</td><td>0.7967</td></tr><tr><td>10</td><td>9π/12</td><td>9π/12</td><td>3π/4</td><td>3π/4</td><td>1.9064</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1.9064/2) = 0.8153</td><td>0.8159</td></tr><tr><td>11</td><td>9π/12</td><td>9π/12</td><td>3π/4</td><td>3π/4</td><td>5.7192</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 3π/4) = 0.9239</td><td>0.9257</td></tr><tr><td>12</td><td>9π/12</td><td>9π/12</td><td>3π/4</td><td>9π/4</td><td>5.7192</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 3π/4) = 0.9239</td><td>0.9257</td></tr><tr><td>13</td><td>9π/12</td><td>9π/12</td><td>9π/4</td><td>9π/4</td><td>5.7192</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 9π/12) = 0.9239</td><td>0.9285</td></tr><tr><td>14</td><td>11π/12</td><td>9π/12</td><td>9π/4</td><td>9π/4</td><td>5.7192</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 9π/12) = 0.9239</td><td>0.9285</td></tr><tr><td>15</td><td>11π/12</td><td>11π/12</td><td>9π/4</td><td>9π/4</td><td>5.7192</td><td>2.7224</td><td>23π/12</td><td>23π/12</td><td>sin(2.7224/2) = 0.9781</td><td>0.9789</td></tr><tr><td>16</td><td>11π/12</td><td>11π/12</td><td>9π/4</td><td>9π/4</td><td>5.7192</td><td>8.1672</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 11π/12) = 0.9914</td><td>0.9946</td></tr><tr><td>17</td><td>11π/12</td><td>13π/12</td><td>9π/4</td><td>9π/4</td><td>5.7192</td><td>8.1672</td><td>23π/12</td><td>23π/12</td><td>sin(1/2 11π/12) = 0.9914</td><td>0.9964</td></tr></tbody></table>
One can see the coefficient $A_k^{(p)}$ as the coefficients in the expansion of a constant function $1/\sqrt{|\partial\Omega|}$ over the orthogonal basis $\{v_k^{(p)}\}$ of $L^2(\partial\Omega)$. As a consequence, one has

$$ \sum_{k=0}^{\infty} |A_k^{(p)}|^2 = 1. \qquad (43) $$

In particular, $|A_k^{(p)}|^2$ can be interpreted as relative weights of different eigenfunctions $v_k^{(p)}$ in some spectral expansions. The convergence of this series implies that $|A_k^{(p)}|$ decreases with $k$, i.e.

$$ \lim_{k \to \infty} A_k^{(p)} = 0. \qquad (44) $$

We inspect the relative contributions of $A_k^{(p)}$ for various domains.

## 4.1 Ellipses and rectangles

In order to break the rotational invariance of a disk, we start again with ellipses of aspect ratio 2 and 10. Figure 5 shows the first 20 coefficients $|A_k^{(p)}|$ for both ellipses. We observe that $A_k^{(p)} \approx 0$ a wide range of $k$ and $p$, except for some coefficients, namely, $A_0^{(p)}$, $A_3^{(p)}$, $A_7^{(p)}$ and $A_{11}^{(p)}$ for the ellipse of aspect ratio 2; and $A_0^{(p)}$, $A_2^{(p)}$ and $A_4^{(p)}$ for the ellipse of aspect ratio 10. Despite the broken rotational invariance of $\Omega$, these results suggest

<figure><img src="image_6.png" /><figcaption>Figure 5: First 20 coefficients $|A_k^{(p)}|$ for (a) an ellipse of semiaxes 1/2 and 1 (aspect ratio 2), (b) an ellipse of semiaxes 1/10 and 1 (aspect ratio 10).</figcaption></figure>

that most of the coefficients $A_k^{(p)}$ vanish. This behavior is a consequence of the domain symmetry. Indeed, as an ellipse has two axes of symmetry (vertical and horizontal), the functions $V_k^{(p)}$ must be symmetric, i.e. to verify $V_k^{(p)}(-x, y) = \pm V_k^{(p)}(x, y)$ and $V_k^{(p)}(x, -y) = \pm V_k^{(p)}(x, y)$. In other words, any $V_k^{(p)}$ satisfies these two relations with one choice of signs from 4 possible combinations. Consequently, the integration of $V_k^{(p)}$ over $\Omega$ yields $A_k^{(p)} = 0$, in 3 out of 4 cases. This is confirmed on Fig. 5a showing the non-zero coefficients $A_0^{(p)}, A_3^{(p)}, A_7^{(p)}$, i.e. 1 out of 4. In the context of diffusion-controlled
reactions, this implies that the related spectral expansions contain only a small number of contributing terms.

After revealing the role of anisotropy of the domain, we inspect to role of the smoothness of the boundary by considering rectangles (Fig. 6). Despite the presence of corners, we still observe that the $A_k^{(p)}$ vanish for most $k$ and a wide range of $p$, except for some coefficients, e.g. $A_0^{(p)}$, $A_5^{(p)}$, and $A_{15}^{(p)}$ for the square, and $A_0^{(p)}$, $A_4^{(p)}$, $A_8^{(p)}$, $A_{10}^{(p)}$ and $A_{16}^{(p)}$ for the considered rectangle. This behavior can be rationalized again by the mirror symmetries of these domains.

<figure><img src="image_7.png" /><figcaption>Figure 6: First 20 coefficients $|A_k^{(p)}|$ for (a) a square of sidelength 2, (b) rectangle with sides 2 and 1.</figcaption></figure>

## 4.2 Prefractals

Figure 7 presents the first 20 coefficients $|A_k^{(p)}|$ for the first three generations of the Koch snowflake. Despite the increased roughness of the prefactal boundary, we still obtain $A_k^{(p)} = 0$ for a wide range of $k$ and $p$, except for $A_0^{(p)}$ and $A_3^{(p)}$ for the equilateral triangle (generation 0); $A_0^{(p)}$, $A_6^{(p)}$ and $A_{18}^{(p)}$ for the generation 1; $A_0^{(p)}$, $A_{12}^{(p)}$ and $A_{18}^{(p)}$ for the generation 2. As previously, cancellation of most coefficients $A_k^{(p)}$ is a consequence of symmetries, which are, however, more sophisticated than in the previous cases of ellipses and rectangles.

## 4.3 Generic triangles

To highlight the role of symmetries in previous examples, we consider a generic triangle with one sidelength 2 and two angles $\pi/3$ and $\pi/12$ shown on Fig. 1d. Figure 7d presents the first 20 coefficients $|A_k^{(p)}|$ but their behavior is very different from the previous ones. Indeed, many coefficients $A_k^{(p)}$ are no longer close to 0. Curiously, the principal eigenfunction $V_0^{(p)}$, which provided the dominant contribution in all previous examples, gives away its leading role to the next eigenfunction $V_1^{(p)}$ for $p \ge 4$.
<figure><img src="image_8.png" /><figcaption>Figure 7: First 20 coefficients $|A_k^{(p)}|$ for (a), (b), (c) the first three generations of the Koch snowflake, and for (d) a triangle with sidelength 2 and angles $\pi/3$ and $\pi/12$.</figcaption></figure>

We conclude that cancellation of many $A_k^{(p)}$ for $p > 0$, which was observed for symmetric domains such as ellipses, rectangles, or even prefractal shapes, is not generic. For an arbitrary domain (e.g. a generic triangle) many terms can contribute to spectral expansions involving $A_k^{(p)}$, even thought their contributions are necessarily reduced as $k \to \infty$. As a consequence, the effect of domain geometry onto various characteristics of diffusion-controlled reactions can be much more sophisticated and versatile than one might expect from earlier theoretical studies focused on highly symmetric domains. Further numerical investigations of this effect and more rigorous characterization of the decay of $|A_k^{(p)}|$ with $k$ present an important perspective for future research.

# 5 Asymptotic behavior of Steklov eigenfunctions

In this section, we analyze the decay of the Steklov eigenfunctions $V_k^{(p)}$ away from the boundary [12, 15, 16, 17, 18]. In particular, Polterovich, Sher and Toth [15] proved that for any bounded domain $\Omega \subseteq \mathbb{R}^2$ with a real-analytic boundary $\partial\Omega$, there exist
positive constants $\eta > 0$ and $B > 0$ depending only on the geometry of $\Omega$, such that all eigenfunctions $V_k^{(0)}$ satisfy

$$ \sqrt{|\partial\Omega|} |V_k^{(0)}(\mathbf{x})| \leq B \exp(-\eta \mu_k^{(0)} |\mathbf{x} - \partial\Omega|), \quad (45) $$

where $|\mathbf{x}-\partial\Omega|$ is the Euclidean distance between $\mathbf{x}$ and the boundary $\partial\Omega$ and we included the prefactor $\sqrt{|\partial\Omega|}$ to ensure that the left-hand side is dimensionless. In other words, each Steklov eigenfunction $V_k^{(0)}$ has an upper bound, which decays exponentially fast away from the boundary, and its decay rate is proportional to the corresponding eigenvalue $\mu_k^{(0)}$. As $\mu_k^{(0)}$ grows with $k$ up to infinity, the exponential decay becomes faster and more and more restrictive. In contrast, the upper bound eq. (45) is in general not much informative for any finite $k$: as $V_k^{(0)}$ is an analytic function on a bounded domain, one can always choose a large enough $B$ or small enough $\eta$ to fulfill the inequality eq. (45). In other words, without restricting estimates on $\eta$ and $B$, the upper bound (45) does not tell much on the behavior of a given Steklov eigenfunction. The upper bound (45) was further generalized to higher dimension and $p \neq 0$ by Helffer and Kachmar [18]. They proved that for any bounded domain $\Omega \subset \mathbb{R}^d$ ($d \geq 2$), with a real-analytic boundary $\partial\Omega$, any eigenfunction $V_k^{(p)}$ with sufficiently high eigenvalue $\mu_k^{(p)}$ (i.e. sufficiently high $k$) has an upper bound that decays exponentially fast away from the boundary. In fact, for any $p > -\lambda_0^D$ (where $\lambda_0^D$ is the smallest eigenvalue of the Dirichlet Laplacian in $\Omega$) there exist constants $B > 0$, $\eta > 0$, $\epsilon > 0$ and $k_0$ such that

$$ \forall k > k_0, \forall \mathbf{x} \in \Omega, \quad |\sqrt{\partial\Omega}||V_k^{(p)}(\mathbf{x})| \leq B \left(\mu_k^{(p)}\right)^{\frac{d}{2}-\frac{1}{4}} \exp\left(-\eta \mu_k^{(p)} \min\{\epsilon, |\mathbf{x} - \partial\Omega|\}\right), \quad (46) $$

where we included again the prefactor $\sqrt{|\partial\Omega|}$. Helffer and Kachmar also questioned whether the assumption of real-analytic boundary could be relaxed.

Being inspired by these fascinating results, we aim to check numerically to which extend the upper bounds eqs. (45) and (46) determine the decay of any given Steklov eigenfunction $V_k^{(p)}$ (with finite $k$). One may wonder whether the exponential decay holds in the whole domain, and whether the eigenvalue $\mu_k^{(p)}$ is the decay rate, i.e. whether $\eta$ is close to 1. To gain intuitive insights onto this behavior, one can first look at the explicit eq. (38) for the Steklov eigenfunctions $V_k^{(p)}$ for the disk of radius $R$. In the limit $p \to 0$ one has for any $k = 1, 2, 3, ...$

$$ \begin{aligned} |V_{2k}^{(0)}(\mathbf{x})| &= \frac{1}{\sqrt{\pi R}} \left(\frac{|\mathbf{x}|}{R}\right)^k = \frac{1}{\sqrt{\pi R}} \exp \left[ k \ln \left( \frac{R - |\mathbf{x} - \partial\Omega|}{R} \right) \right] \\ &\leq \frac{1}{\sqrt{\pi R}} \exp \left( -\frac{k}{R} |\mathbf{x} - \partial\Omega| \right) = \frac{\exp \left( -\mu_{2k}^{(0)} |\mathbf{x} - \partial\Omega| \right)}{\sqrt{\pi R}}, \quad (47) \end{aligned} $$

where we used that $\mu_{2k}^{(0)} = k/R$; and the same applies for $V_{2k+1}^{(0)}(\mathbf{x})$. One sees how $\mu_k^{(0)}$ controls the exponential decay of $V_k^{(0)}$ away from the boundary, with $\eta = 1$ and $B = 1/\sqrt{\pi R}$. In the opposite limit $p \to \infty$, the asymptotic behavior of the modified
Bessel functions $I_k(z)$ implies for $r\sqrt{p} \gg 1$:

$$ |V_k^{(p)}(\mathbf{x})| \lesssim \frac{1}{\sqrt{\pi R}} \frac{\sqrt{R}}{\sqrt{r}} \exp(-\sqrt{p}(R-r)) \approx \frac{1}{\sqrt{\pi R}} \frac{\sqrt{R}}{\sqrt{r}} \exp(-\mu_k^{(p)}|\mathbf{x} - \partial\Omega|) \quad (48) \\ \leq \frac{p^{1/4}}{\sqrt{\pi}} \exp(-\mu_k^{(p)}|\mathbf{x} - \partial\Omega|), $$

where we used the asymptotic relation (9). In turn, if $r$ is so small that $r\sqrt{p} \ll 1$ (but still $R\sqrt{p} \gg 1$) one uses $I_k(r\sqrt{p}) \approx \frac{(r\sqrt{p})^k}{k!}$ to get

$$ |V_{2k}^{(p)}(\mathbf{x})| \ll \frac{p^{1/4}}{k!} \exp(-\sqrt{p}R) \lesssim \frac{p^{1/4}}{k!} \exp(-\mu_{2k}^{(p)}|\mathbf{x} - \partial\Omega|). \quad (49) $$

In both limits, we observe an exponential decay of $V_k^{(p)}$ away from the boundary, which is controlled by $\mu_k^{(p)}$, and $\eta$ is close to 1 (even though the above upper bounds suggest $\eta = 1$, our derivation involved some asymptotic relations, which may require having $\eta$ slightly smaller than 1). In Appendix 7.2, we analyze the explicit representation of the Steklov eigenfunctions for a rectangle and come to the same conclusion for $p = 0$. In sharp contrast to the disk, the boundary of the rectangle is not real-analytic due to the corners. Nevertheless, the upper bound eq. (45) holds and the eigenfunctions exhibit exponential decay controlled by $\mu_k^{(0)}$, with $\eta$ close to 1.

Let us give another argument in favor of the particular choice $\eta = 1$. Let $\mathbf{x}_0$ be a boundary point, and $\mathbf{x}_\delta = \mathbf{x}_0 - \delta n_{\mathbf{x}_0} \in \Omega$ be a bulk point at small distance $\delta$ from $\partial\Omega$, where $n_{\mathbf{x}_0}$ is the normal unit vector to $\partial\Omega$ at $\mathbf{x}_0$ oriented outward the domain. If we assume that $V_k^{(p)}(\mathbf{x}_0) > 0$, then the Steklov boundary condition can be written as

$$ \partial_n \ln (V_k^{(p)}(\mathbf{x}_0)) = \mu_k^{(p)} \quad (\mathbf{x}_0 \in \partial\Omega) \quad (50) $$

(if $V_k^{(p)}(\mathbf{x}_0) < 0$, one can replace $V_k^{(p)}(\mathbf{x}_0)$ by $-V_k^{(p)}(\mathbf{x}_0)$ that yields $\partial_n \ln (|V_k^{(p)}(\mathbf{x}_0)|)|_{\partial\Omega} = \mu_k^{(p)}$, with no changes in the argument; in turn, we ignore here the specific points at which $V_k^{(p)}(\mathbf{x}_0) = 0$). Since $V_k^{(p)}$ is analytic in $\Omega$ and $V_k^{(p)}(\mathbf{x}_0) > 0$, it is also positive in a small vicinity of the boundary, so that

$$ \ln (V_k^{(p)}(\mathbf{x}_\delta)) = \ln (V_k^{(p)}(\mathbf{x}_0)) - \mu_k^{(p)}\delta + O(\delta^2), \quad (51) $$

to be consistent with eq. (50). As a consequence, the Steklov eigenfunction $V_k^{(p)}$ exhibits an exponential decay near the boundary

$$ V_k^{(p)}(\mathbf{x}_\delta) \approx V_k^{(p)}(\mathbf{x}_0) \exp(-\mu_k^{(p)}\delta), \quad (52) $$

where $\delta = |\mathbf{x}_\delta - \partial\Omega|$ is the distance to the boundary. One sees that the value $\eta = 1$ naturally comes from the Steklov condition. The fundamental question is whether the exponential behavior (52) holds approximately far from the boundary, when the above argument is not applicable.
We address this question numerically for $p = 0$ and check the exponential decay of $V_k^{(0)}$ away from the boundary for various planar domains, even when $\partial\Omega$ is not real-analytic. Most importantly, we discuss whether $\eta$ is close to 1 in general. The top row of Fig. 8 presents the log-scaled eigenfunction $V_{15}^{(0)}$ for three polygonal domains. It confirms the expected exponential decay away from the boundary for the square (see Appendix 7.2), but also shows the exponential decay of $V_{15}^{(0)}$ for a pentagon and the second generation of the Koch snowflake. In these cases, the presence of corners does not seem to affect the localization near the smooth parts of the polygonal boundary. For a more systematic insight, we define

$$B_k(\mathbf{x}) = \left| \sqrt{|\partial\Omega|} V_k^{(0)}(\mathbf{x}) \exp\left(\mu_k^{(0)}|\mathbf{x} - \partial\Omega|\right) \right|, \qquad (53)$$

as if we explicitly set $\eta = 1$ in eq. (45). If $B_k(\mathbf{x})$ was constant, the eigenfunction $V_k^{(0)}(\mathbf{x})$ would exhibit the exponential decay away from the boundary with the rate $\mu_k^{(0)}$. As a consequence, variations of $B_k(\mathbf{x})$ and, in particular, its high values can indicate regions where the exponential decay $\exp(-\mu_k^{(0)}|\mathbf{x} - \partial\Omega|)$ does not hold. We note that, even for a disk, oscillations of $V_k^{(0)}(\mathbf{x})$ on the boundary and inside the domain result in variations of $B_k(\mathbf{x})$. The bottom row of Fig. 8 presents the function $B_{15}(\mathbf{x})$ for the three polygonal domains. For instance, we get $\max_{\mathbf{x}\in\Omega} B_{15}(\mathbf{x}) \approx 5.6$, $\max_{\mathbf{x}\in\Omega} B_{15}(\mathbf{x}) \approx 2.6$ and $\max_{\mathbf{x}\in\Omega} B_{15}(\mathbf{x}) \approx 3.8$ for the square, the pentagon and the Koch snowflake respectively. Interestingly, the maximal deviation from the upper bound can be found in the center (for a square), in the central part (for a pentagon), or near the boundary (for a Koch snowflake). However, for the pentagon and the Koch snowflake, we also observe that there exist indices $k$ for which the maximum value of $B_k(\mathbf{x})$ becomes much larger, suggesting deviations of the related eigenfunctions from the exponential decay with $\eta$ close to 1.

To better understand the origin of this behavior, we go back to real-analytic boundaries and focus on the disk and its minor deformations. Using polar coordinates $(r, \theta)$, we define the boundary of the deformed disks $\Omega_\gamma = \{(r, \theta) | r < \rho(\theta)\}$ by setting

$$\rho(\theta) = 1 + \gamma \cos(5\theta), \quad 0 \le \theta \le 2\pi. \qquad (54)$$

Figure 9 presents the localization and decay of $V_{20}^{(0)}$ away from the boundary, as well as the function $B_{20}(\mathbf{x})$ for a disk, and the two slightly deformed disks with $\gamma = 0.01$ and $\gamma = 0.02$. As already observed in [33], a very small perturbation of the disk may result in significant changes of the structure of nodal lines and thus in the interior behavior of Steklov eigenfunctions. Indeed, a very small perturbation results in high values of $B_k(\mathbf{x})$ at the center of the domain for some specific indices $k$, e.g. for $k = 20$. These high values of $B_k(\mathbf{x})$ suggest that the related eigenfunction decreases slower inside the domain than near the boundary; in other words, the upper bound (45) with $\eta \approx 1$ may not be valid for the whole domain. Curiously, we do not retrieve such behavior for ellipses (see Fig. 10), which could also be seen as deformations of a disk. This can be a consequence of symmetries.

In order to investigate the exponential decay of eigenfunctions in the whole domain, we introduce the function

$$U_k^{(p)}(\delta) = \sqrt{|\partial\Omega|} \max_{\mathbf{x}\in\gamma_\delta} |V_k^{(p)}(\mathbf{x})|, \qquad (55)$$
<figure><img src="image_9.png" /><figcaption>Figure 8: (Top row) Localization and decay of $V_{15}^{(0)}$ away from the boundary for (a) a square, (b) a pentagon, (c) generation 2 of the Koch snowflake. Colormaps illustrate the behavior of $\lg(|V_{15}^{(0)}|)$, where $\lg(z) = \ln(z)/\ln(10)$ is the decimal logarithm; white regions correspond to the values $|V_{15}^{(0)}| < 10^{-4}$. The related eigenvalues are (a) $\mu_{15}^{(0)} \approx 5.50$, (b) $\mu_{15}^{(0)} \approx 9.00$ and (c) $\mu_{15}^{(0)} \approx 2.75$. (Bottom row) The function $B_{15}(x)$ from eq. (53).</figcaption></figure>

where $\gamma_{\delta}$ is the contour line of points in $\Omega$ at distance $\delta$ from the boundary $\partial\Omega$. We aim to test whether the following approximation holds:

$$U_k^{(p)}(\delta) \approx U_k^{(p)}(0)e^{-\mu_k^{(p)}\delta}. \qquad (56)$$

Figure 11a shows the log-plot of $U_{20}^{(0)}(\delta)$ for the three eigenfunctions shown on Fig. 9 and compare its exponential decay with that known for the disk. One sees that the more the disk is perturbed, the earlier the exponential decay $\exp(-\mu_k^{(0)}\delta)$ stops to approximate the behavior of the related eigenfunction $V_k^{(0)}$. This indicates that the determination of an optimal $\eta$ in eq. (45) might not be straightforward even for domains with real-analytic boundary and requires further analysis. Note that the numerical results shown in Fig. 11 are more accurately described by the truncated exponential $\exp(-\mu_k^{(0)} \min\{\epsilon, |x - \partial\Omega|\})$ that appears in eq. (46). Indeed, there is a cutoff distance $\epsilon$ above which the exponential decay with the decay rate $\mu_k^{(0)}$ does not work. In other words, $\epsilon$ determines an inner “central” region of the domain in which the eigenfunction decays slower. In turn, Fig. 11b illustrates that there is no such a central region for the considered ellipses, i.e., there may not be need in the truncated exponential for these domains. From a practical point of view, the natural question is to know whether the exponential function
<figure><img src="image_10.png" /><figcaption>Figure 9: (Top row) Localization and decay of $V_{20}^{(0)}$ away from the boundary for (a) a disk, and two slightly deformed disks defined by eq. (54) with (b) $\gamma = 0.01$, and (c) $\gamma = 0.02$. Colormaps illustrate the behavior of $\lg(|V_{20}^{(0)}|)$, where $\lg(z) = \ln(z)/\ln(10)$ is the decimal logarithm; white regions correspond to the values $|V_{20}^{(0)}| < 10^{-4}$. The related eigenvalues are (a) $\mu_{20}^{(0)} = 10.00$, (b) $\mu_{20}^{(0)} \approx 9.99$ and (c) $\mu_{20}^{(0)} \approx 9.98$. (Bottom row) The function $B_{20}(\mathbf{x})$ from eq. (53).</figcaption></figure>

$b_k \exp(-\mu_k^{(0)}|\mathbf{x} - \partial\Omega|)$, or its truncated form $b_k \exp(-\mu_k^{(0)} \min\{\epsilon, |\mathbf{x} - \partial\Omega|\})$, can be an accurate approximation of the Steklov eigenfunction $V_k^{(0)}$. For instance, can one choose $b_k$ such that $\max_{\mathbf{x} \in \Omega} |V_k^{(p)}(\mathbf{x}) - b_k \exp(-\mu_k^{(0)}|\mathbf{x} - \partial\Omega|)|$ is small enough?

Finally, we complete this section by providing a complementary insight onto the localization of $V_k^{(p)}$ at large $p$. Multiplying eq. (7) by $V_k^{(p)}$, integrating over $\Omega$, using the Green's formula, and employing the Steklov boundary condition, one easily gets

$$p \int_{\Omega} [V_k^{(p)}(\mathbf{x})]^2 d\mathbf{x} + \int_{\Omega} |\nabla V_k^{(p)}(\mathbf{x})|^2 d\mathbf{x} = \mu_k^{(p)}. \qquad (57)$$

The eigenvalue determines therefore the combination of $L^2(\Omega)$-norms of $V_k^{(p)}$ and of its gradient. It is instructive to compute both norms separately. In Appendix 7.3, we derive the following relation

$$\int_{\Omega} |V_k^{(p)}(\mathbf{x})|^2 d\mathbf{x} = \partial_p \mu_k^{(p)}. \qquad (58)$$
<figure><img src="image_11.png" /><figcaption>Figure 10: (Top row) Localization and decay of $V_{20}^{(0)}$ away from the boundary for (a) an ellipse with semiaxes $a = 1$, $b = 1.01$ (b) an ellipse with semiaxes $a = 1$, $b = 0.8$, and (c) an ellipse with semiaxes $a = 1$, $b = 0.5$. Colormaps illustrate the behavior of $\lg(|V_{20}^{(0)}|)$. (Bottom row) The function $B_{20}(\boldsymbol{x})$ from eq. (53). The related eigenvalues are (a) $\mu_{20}^{(0)} \approx 9.95$, (b) $\mu_{20}^{(0)} \approx 11.08$ and (c) $\mu_{20}^{(0)} \approx 12.98$.</figcaption></figure>

Combining both equations, we also find

$$ \int_{\Omega} |\nabla V_k^{(p)}(\boldsymbol{x})|^2 d\boldsymbol{x} = \mu_k^{(p)} - p \partial_p \mu_k^{(p)}. \qquad (59) $$

As we discussed in Sec. 3, the eigenvalues $\mu_k^{(p)}$ grow as $\sqrt{p}$ at large $p$, and so does the left-hand side of Eq. (57). In turn, Eq. (58) implies that the $L^2(\Omega)$-norm of $V_k^{(p)}$ vanishes as $p^{-1/4}$ as $p \to \infty$. This is a weaker form of the localization of the Steklov eigenfunctions near the boundary.

# 6 Discussion and conclusion

In this paper, we numerically investigated the spectral properties of the Dirichlet-to-Neumann operator $M_p$ and their dependence on the parameter $p$ and the domain geometry. We considered various shapes, including ellipses, triangles, rectangles, regular polygons and prefractal Koch snowflakes.

Our first contribution concerned the asymptotic behavior of the eigenvalues of $M_p$. For all considered shapes, we confirmed the validity of the asymptotic relations eqs. (10) and (40) for large and small $p$. In the limit $p \to \infty$, the coefficients $c_k$ in eq. (10) were
<figure><img src="image_12.png" /><figcaption>Figure 11: Log-plot of the decay of $U_{20}^{(0)}(\delta)$ for (a) a disk ($\gamma = 0$) and two deformed disks ($\gamma = 0.01$ and $\gamma = 0.02$), where solid black line indicates the known decay $\sqrt{2}\exp(-10\delta)$ for a disk; (b) an ellipse with semiaxes $a = 1$, $b = 1.01$, an ellipse with semiaxes $a = 1$, $b = 0.8$, and an ellipse with semiaxes $a = 1$, $b = 0.5$, where lines indicate the expected decay $U_{20}^{(0)}(0)\exp(-\mu_{20}^{(0)}\delta)$. The related eigenvalues are $\mu_{20}^{(0)} \approx 9.95$, $\mu_{20}^{(0)} \approx 11.08$ and $\mu_{20}^{(0)} \approx 12.98$ respectively.</figcaption></figure>

known to be 1 for bounded domains with $C^1$ boundary. In turn, their values for polygonal domains were unknown. We conjectured that as $p$ increases, first Steklov eigenfunctions are getting localized near the corners of a polygonal domain, and suggested an iterative procedure to obtain the coefficients $c_k$. This conjectural relation was numerically validated on several domains. We are unaware of earlier studies on the coefficients $c_k$ for the asymptotic behavior of the eigenvalues $\mu_k^{(p)}$ of the Dirichlet-to-Neumann operator. However, additional insights can be gained from the related problem of the Robin Laplacian defined in eq. (36). In fact, Lacey et al. considered the asymptotic behavior of the smallest eigenvalue $\lambda_0^{(q)}$ of the Robin Laplacian in the limit $q \to -\infty$, and found that [39]:

$$ \lambda_0^{(q)} \approx -q^2 / \sin^2(\alpha/2) \qquad (60) $$

in a corner with angle $\alpha < \pi$. This result was generalized by Levitin and Parnovski [40] who proved eq. (60) for a polygonal domain with angles $\alpha_1, \dots, \alpha_n$, such that $0 < \alpha_i < \pi$, and $\alpha = \min\{\alpha_i\}$. Then, Khalile [41] extended this result for the first $n$ eigenvalues of the Robin Laplacian for polygons with angles between $\pi/3$ and $\pi$. Finally, Khalile and Pankrashkin [42] considered a similar problem for an infinite sector and showed the asymptotic behavior $\lambda_n^{(q)} \approx -q^2/((2n+1)\alpha/2)^2$ in the limit $\alpha \to 0$. The duality between the Robin problem and the Steklov problem allows one to invert eq. (60) to get $\mu_0^{(p)} \approx \sin(\alpha/2)\sqrt{p}$ as $p \to \infty$ for a polygon. This asymptotic result, which is a direct consequence of [40], confirms our conjecture for the smallest eigenvalue $\mu_0^{(p)}$. In turn, our conjectural iterative procedure relates the asymptotic behavior of all eigenvalues to the angles of a polygonal domain. Its rigorous demonstration presents an interesting open problem. Perhaps, the most challenging part is to prove that our procedure yields all coefficients $c_k < 1$.
The second result concerned some spectral expansions that appear in the theory of diffusion-controlled reactions [10]. For symmetric domains such as ellipses, rectangles and even prefactual Koch snowflakes, we numerically observed that many coefficients $A_k^{(p)}$ vanish for a wide range of $k$ and $p$. We argued that this behavior is a consequence of domain symmetries. In turn, this property breaks for generic domains such as an arbitrary triangle. On one hand, the cancellation of many coefficients $A_k^{(p)}$ due to domain symmetries can considerably simplify spectral expansions and the analysis of diffusion-controlled reactions in such domains. Moreover, these symmetries can potentially be used to design domains with specific properties. On the other hand, the use of symmetric domains as examples in theoretical and numerical studies may lead to erroneous conjectures and conclusions. For instance, the cancellation of coefficients $A_k^{(p)}$ with $k > 0$ for the disk is a very specific consequence of its rotational symmetry, which fails for less symmetric domains. In particular, many coefficients $A_k^{(p)}$ contribute to the spectral expansion in the case of a generic triangle. It is therefore an open question how many terms are relevant and how their number depends on domain shape. More formally, it is known from the general arguments that $A_k^{(p)} \to 0$ as $k \to \infty$, but the speed of decay remains unknown. This question is practically important because a rapid decay of $A_k^{(p)}$ may allow truncating some spectral expansions to get useful approximations with a limited number of contributing eigenmodes.

Last, we investigated the localization of Steklov eigenfunctions $V_k^{(p)}$ on the boundary in the presence of corners. Even if the boundary is polygonal and thus is not real-analytic, Fig. 8 illustrated the exponential decay of $|V_k^{(0)}|$ away from the boundary for a broad range of indices $k$. These numerical examples could motivate further mathematical investigations in this direction. However, we also observed that there exist domains and indices $k$, for which deviations from the exponential decay, characterized by $B_k(x)$ and $U_k^{(p)}(\delta)$, become large, suggesting that the hypothesis $\eta \approx 1$ in eq. (45) might not be relevant for all $k$ and in the whole domain. In particular, we noticed that only a slight perturbation of the disk can result in high values of $B_k(x)$ located in the center of the domain, indicating that the related eigenfunctions decrease slower inside the domain than near the boundary. A more systematic study of this behavior through the estimation of the constants $B$ and $\eta$ in eq. (45) and their relation with the domain geometry present an interesting perspective to this work.

## Acknowledgments

The authors thank Prof. Iosif Polterovich, Prof. Michael Levitin and Prof. Bernard Hellffer for fruitful discussions and suggestions. The authors acknowledge Antoine Moutal for his early contributions in the numerical implementation of the finite element method.

# 7 Appendix

## 7.1 Numerical validation: solution for rectangles

Laugesen studied the spectral properties of the Robin Laplacian in various rectangular domains (or cuboids) in $\mathbb{R}^d$ [43]. In this section, we adapt his analysis to the Steklov
problem and give the explicit formulas for rectangles $\Omega = (0, b_1) \times (0, b_2) \subset \mathbb{R}^2$. The separation of variables in the modified Helmholtz equation $(p - \Delta)u = 0$ yields:

$$u(\boldsymbol{x}) = \prod_{n=1}^{2} \left( \alpha_n \cosh(\alpha_n x_n / b_n) - \mu b_n \sinh(\alpha_n x_n / b_n) \right), \qquad (61)$$

with

$$\tanh(\alpha_1) = \frac{2\mu\alpha_1 b_1}{\alpha_1^2 + \mu^2 b_1^2}, \qquad (62)$$

$$\tanh(\alpha_2) = \frac{2\mu\alpha_2 b_2}{\alpha_2^2 + \mu^2 b_2^2}, \qquad (63)$$

and

$$\frac{\alpha_1^2}{b_1^2} + \frac{\alpha_2^2}{b_2^2} = p. \qquad (64)$$

Solving the system of three nonlinear equations (62), (63), (64), one can determine the unknown coefficients $\alpha_1, \alpha_2$ and the eigenvalue $\mu$. We stress that solutions $\alpha_n$ can be either real, or purely imaginary. In order to solve the system of three nonlinear equations, we consider eq. (62) as the quadratic equation on $\mu$, whose two solutions are

$$\mu_+ = \frac{\alpha_1}{b_1} \tanh(\alpha_1/2), \quad \mu_- = \frac{\alpha_1}{b_1} \tanh(\alpha_1/2). \qquad (65)$$

Let us first focus on the case $p = 0$, for which eq. (64) implies $\alpha_2 = i\alpha_1 b_2 / b_1$. Substituting this expression into eq. (63), we get

$$\tan(\alpha_1 b_2 / b_1) = \frac{2\alpha_1 b_1 \mu}{\mu^2 b_1^2 - \alpha_1^2}. \qquad (66)$$

Substituting $\mu_{\pm}$ from Eq. (65) into this relation, one has

$$\tan(\alpha_1 b_2 / b_1) = \pm \sinh(\alpha_1). \qquad (67)$$

Each of these equations has infinitely many real solutions that determine the eigenvalues of the Dirichlet-to-Neumann operator $M_0$ for a rectangle.

When $b_1 = b_2 = 2$, one can set $\kappa = \alpha_1/2$ and use trigonometric relations to rewrite Eq. (67) with plus sign in terms of $\kappa$ as:

$$\frac{\tan(\kappa)}{1 - \tan^2(\kappa)} = \frac{\tanh(\kappa)}{1 - \tanh^2(\kappa)} = \frac{-\tanh(\kappa)}{1 - \tanh^2(\kappa)},$$

which is equivalent to $\tan(\kappa) = \tanh(\kappa)$ or $\tan(\kappa) = -\tanh(\kappa)$. Similarly, eq. (67) with minus sign reads

$$\frac{\tan(\kappa)}{1 - \tan^2(\kappa)} = \frac{\tanh(\kappa)}{1 - \tanh^2(\kappa)} = \frac{-\tanh(\kappa)}{1 - \tanh^2(\kappa)},$$

which is equivalent to $\tan(\kappa) = \tanh(\kappa)$ or $\tan(\kappa) = -\tanh(\kappa)$. We retrieve therefore the equations reported in Table 7.1 of [4] for the square $(-1, 1)^2$. In this case, there is an
additional eigenvalue 1, which corresponds to the eigenfunction $x_1x_2$. This eigenfunction is not included in the general form (61).

Now we return to the case $p > 0$. We search for real solutions $\alpha_1$. Setting $\alpha_1 = (b_1/b_2)\sqrt{-\alpha_2^2 + pb_2^2}$ and substituting $\mu_+$ from Eq. (65), we get the following equation on $\alpha_2$:

$$ \tanh(\alpha_2) = -\alpha_2 \frac{(\alpha_1 b_2 / b_1) \sinh(\alpha_1)}{\alpha_2^2 - pb_2^2 \cosh^2(\alpha_1/2)}. \qquad (68) $$

There is a finite number of real solutions $\alpha_2$ of this equation on the interval from 0 to $b_2\sqrt{-\lambda}$, for which $\alpha_1$ is real. In turn, there are infinitely many purely imaginary solutions $\alpha_2$. Setting $\alpha_2 = i\alpha$ and thus $\alpha_1 = (b_1/b_2)\sqrt{\alpha_2^2 + pb_2^2}$, we transform the above equation into

$$ \tan(\alpha) = \alpha \frac{(\alpha_1 b_2 / b_1) \sinh(\alpha_1)}{\alpha^2 + pb_2^2 \cosh^2(\alpha_1/2)}. \qquad (69) $$

As $p > 0$, $\alpha_1$ is real and the right-hand side of Eq. (69) is positive. This equation has infinitely many solutions, which lie on the intervals $(\pi k, \pi k + \pi/2)$, with $k = 0, 1, 2, \dots$.

Similarly, using $\mu_-$ from Eq. (65), we get another equation on $\alpha_2$:

$$ \tanh(\alpha_2) = \alpha_2 \frac{(\alpha_1 b_2 / b_1) \sinh(\alpha_1)}{\alpha^2 + pb_2^2 \sinh^2(\alpha_1/2)}, \qquad (70) $$

which has a finite number of real solutions $\alpha_2$ on the interval $(0, b_2\sqrt{p})$. In turn, setting again $\alpha_2 = i\alpha$, one gets

$$ \tan(\alpha) = -\alpha \frac{(\alpha_1 b_2 / b_1) \sinh(\alpha_1)}{\alpha^2 - pb_2^2 \sinh^2(\alpha_1/2)}. $$

Since the denominator in the right-hand side can change sign, it is more convenient to rewrite this relation as

$$ \frac{\alpha}{\tan(\alpha)} = - \frac{\alpha^2 - pb_2^2 \sinh^2(\alpha_1/2)}{(\alpha_1 b_2 / b_1) \sinh(\alpha_1)}. \qquad (71) $$

There are infinitely many solutions of this equation.

The above computation allows one to find all solutions, for which $\alpha_1$ is real, while $\alpha_2$ is either real, or purely imaginary. Exchanging the roles of $\alpha_1$ and $\alpha_2$, one can also determine the missing pairs, for which $\alpha_2$ is real while $\alpha_1$ is either real, or purely imaginary. Combining all these solutions, we determine the eigenvalues of the Dirichlet-to-Neumann operator $M_p$ for the rectangle.

Table 5 summarizes the first 11 eigenvalues of $M_p$ for a rectangles of sides 1 and 2, with $p = 1$ (Fig. 1e shows the associated eigenfunction $V_4^{(1)}$). One sees that the eigenvalues in the third column, which were numerically obtained by Method 1, are in excellent agreement with the exact ones.

Let us inspect the limit $p \to \infty$. Setting $0 \le \gamma \le \pi/2$, it is convenient to write $\alpha_1 = b_1\sqrt{p} \sin \gamma$ and $\alpha_2 = b_2\sqrt{p} \cos \gamma$ that satisfy Eq. (64). Substituting these expressions into Eq. (68), we get

$$ \tanh(b_2\sqrt{p} \cos \gamma) = \frac{\sin(2\gamma) \sinh(b_1\sqrt{p} \sin \gamma)}{2[\cosh^2(b_1\sqrt{p} \sin \gamma/2) - \cos^2(\gamma)]}. \qquad (72) $$
<table><caption>Table 5: List of the first 11 eigenvalues $\mu_k^{(p)}$ for a rectangle with sides 1 and 2, with $p = 1$. For Method 1, the mesh is composed of 179918 triangles, and the maximal mesh size is 0.005. For Method 2, the mesh is composed of 16256 triangles, the maximum mesh size is 0.03, and the series in eq. (35) was truncated to 88 eigenfunctions of the Laplace operator.</caption><thead><tr><th>Index k</th><th>Exact</th><th>Method 1</th><th>Method 2</th></tr></thead><tbody><tr><td>0</td><td>0.3105</td><td>0.3105</td><td>0.3105</td></tr><tr><td>1</td><td>0.7511</td><td>0.7511</td><td>0.7512</td></tr><tr><td>2</td><td>1.6451</td><td>1.6451</td><td>1.6435</td></tr><tr><td>3</td><td>1.7342</td><td>1.7342</td><td>1.7332</td></tr><tr><td>4</td><td>2.2304</td><td>2.2305</td><td>2.2315</td></tr><tr><td>5</td><td>2.9051</td><td>2.9051</td><td>2.9219</td></tr><tr><td>6</td><td>3.9156</td><td>3.9159</td><td>3.9232</td></tr><tr><td>7</td><td>4.1665</td><td>4.1668</td><td>4.2058</td></tr><tr><td>8</td><td>4.7951</td><td>4.7955</td><td>4.8040</td></tr><tr><td>9</td><td>4.7961</td><td>4.7966</td><td>4.8147</td></tr><tr><td>10</td><td>5.5419</td><td>5.5426</td><td>5.5582</td></tr></tbody></table>

One can easily check that $\gamma \approx 0$ is not compatible with this equation. As a consequence, since $b_1\sqrt{p}\sin\gamma \gg 1$ in the limit $p \to \infty$, the right-hand side of the above equation is close to $\sin(2\gamma)$. If $b_2\sqrt{p}\cos\gamma \gg 1$, one gets $1 \approx \sin(2\gamma)$, from which $\gamma = \pi/4$ and thus $\mu \approx \alpha_1/b_1 \approx \sin(\pi/4)\sqrt{p}$. In turn, if $b_2\sqrt{p}\cos\gamma \ll 1$, then $\gamma \approx \pi/2$, and thus $\mu \approx \sqrt{p}$. These qualitative arguments are consistent with our numerical predictions on the behavior of the eigenvalues in the limit $p \to \infty$. However, more accurate analysis is needed to claim that only two eigenvalues behave as $\mu \approx \sin(\pi/4)\sqrt{p}$.

## 7.2 Localization in rectangles

Let us focus on $p = 0$ and consider the solutions of eq. (66) for which $\alpha_1$ is real, while $\alpha_2 = i\alpha_1 b_2/b_1$ is purely imaginary. The associated eigenfunction, which is given by eq. (61) up to a normalization, is factored as $u_1(x_1)u_2(x_2)$. As the factor $u_2(x_2)$ exhibits oscillatory behavior along $x_2$ coordinate, we focus on $u_1(x_1)$, which can be written as

$$u_1(x_1) = \frac{\alpha_1 - \mu b_1}{2} e^{\alpha_1 x_1/b_1} + \frac{\alpha_1 + \mu b_1}{2} e^{-\alpha_1 x_1/b_1}. \quad (73)$$

Using the first relation in eq. (65) to express $\mu$, one gets then

$$u_1(x_1) = \frac{\alpha_1}{1 - e^{-\alpha_1}} \left( e^{-\alpha_1 x_1/b_1} - e^{-\alpha_1(1-x_1/b_1)} \right), \quad (74)$$

i.e., $|u_1(x_1)|$ decays exponentially with the distance $\delta_1 = \min\{x_1, b_1 - x_1\}$ from either of two endpoints of the interval $(0, b_1)$:

$$|u_1(x_1)| \le \frac{\alpha_1}{1 - e^{-\alpha_1}} e^{-\delta_1 \alpha_1 / b_1} \approx \frac{\alpha_1}{1 - e^{-\alpha_1}} e^{-\mu_+ \delta_1}, \quad (75)$$

where $\alpha_1/b_1$ is exponentially close to $\mu_+$ when $\alpha_1$ is large enough.
Similarly, if one uses the second relation in eq. (65) to express $\mu$, one gets

$$u_1(x_1) = \frac{\alpha_1}{1 + e^{-\alpha_1}} \left( e^{-\alpha_1 x_1 / b_1} + e^{-\alpha_1 (1 - x_1 / b_1)} \right), \quad (76)$$

which can be bounded as

$$|u_1(x_1)| \le \frac{2\alpha_1}{1 + e^{-\alpha_1}} e^{-\delta_1 \alpha_1 / b_1} \approx \frac{2\alpha_1}{1 + e^{-\alpha_1}} e^{-\mu - \delta_1}. \quad (77)$$

Setting $\delta_2 = \min\{x_2, b_2 - x_2\}$, one gets $|x - \partial\Omega| = \min\{\delta_1, \delta_2\} \le \delta_1$, so that the upper bounds eqs. (75) and (77) can be written as

$$|u_1(x_1)| \le C_1 e^{-\mu |x - \partial\Omega|}, \quad (78)$$

with a constant $C_1$. Finally, the analytic function $u_2(x_2)$ is bounded by its maximum, one has

$$|\sqrt{\partial\Omega}||V_k^{(0)}| \lesssim B \exp(-\mu_k^{(0)}|x - \partial\Omega|), \quad (79)$$

i.e. we retrieve the upper bound (45) with $\eta \approx 1$.

We note that eqs. (74) and (76) highlight the expected symmetry of the eigenfunctions discussed in Sec. 4, namely, $u_1(x_1) = -u_1(b_1 - x_1)$ in eq. (74) and $u_1(x_1) = u_1(b_1 - x_1)$ in eq. (76). In the former case, the antisymmetric eigenfunction vanishes in the middle so that it may decay even faster than exponential near the middle point. In turn, the symmetric eigenfunction does not vanish at the middle, and both terms in eq. (76) provide equal contributions in the middle. Such a function is expected to decay slower in the middle. Figure 12a illustrates this behavior for eigenfunctions $V_{10}^{(0)}$ and $V_{11}^{(0)}$. One sees that, even though the exponentially decaying upper bound eq. (79) holds for both cases, it does not necessarily approximate the eigenfunction. In fact, to fulfill the upper bound for symmetric eigenfunctions, we had to add the factor 2 in eq. (77), which shifts the upper bound from the expected exponential behavior of the eigenfunction. The above symmetry argument is applicable to eigenfunctions corresponding to simple (non-degenerate) eigenvalues. In turn, if two (or more) eigenfunctions correspond to the same eigenvalue, then their linear superposition is also an eigenfunction that can break this symmetry. We illustrate this situation for a rectangle with sides 2 and 1, for which Figs. 12b and 12c shows two eigenfunctions $V_{20}^{(0)}$ and $V_{21}^{(0)}$ that correspond to a twice degenerate eigenvalue. One of these eigenfunctions is localized on the left edge and the other is on the right edge of the rectangle.

Refining the above arguments, one can achieve more rigorous statements, in particular on the closeness of $\eta$ to 1. We expect that similar analysis can be performed for $p > 0$.

## 7.3 Derivation of the $L^2(\Omega)$-norm of $V_k^{(p)}$

In this Appendix, we derive the identity (58). First, applying the spectral expansion eq. (35) and the orthogonality of Laplacian eigenfunctions $u_k^{(q)}$ to each other, we get a simple identity for Green's functions

$$\int_{\Omega} \tilde{G}_q(\mathbf{x}_1, p|\mathbf{x}) \tilde{G}_q(\mathbf{x}_2, p|\mathbf{x}) d\mathbf{x} = -\partial_p \tilde{G}_q(\mathbf{x}_2, p|\mathbf{x}_1) \quad (\forall \mathbf{x}_1, \mathbf{x}_2 \in \Omega). \quad (80)$$
<figure><img src="image_13.png" /><figcaption>Figure 12: (a) Log-plot of $U_{10}^{(0)}(\delta)$ and $U_{11}^{(0)}(\delta)$ for a rectangle with sides 2 and 1, the related eigenvalues are $\mu_{10}^{(0)} \approx 5.46$ and $\mu_{11}^{(0)} \approx 5.54$, and lines indicate the expected decay $U_{10}^{(0)}(0) \exp(-\mu_{10}^{(0)}\delta)$ and the upper bound $2 \cdot U_{10}^{(0)}(0) \exp(-\mu_{10}^{(0)}\delta)$ for one of the eigenfunctions; (b) and (c) the localization and decay of $V_{20}^{(0)}$ and $V_{21}^{(0)}$ away from the boundary. Colormaps illustrate the behavior of $\lg(|V_{20}^{(0)}|)$ on (b), and $\lg(|V_{21}^{(0)}|)$ on (c). The related eigenvalues are $\mu_{20}^{(0)} = \mu_{21}^{(0)} \approx 11.00$.</figcaption></figure>

Next, using the representation (37), one can write the squared $L^2(\Omega)$-norm of the Steklov eigenfunction $V_k^{(p)}$ as

$$ \begin{aligned} \|V_k^{(p)}\|_{L^2(\Omega)}^2 &= \int_{\Omega} |V_k^{(p)}(\mathbf{x})|^2 d\mathbf{x} = \int_{\partial\Omega} \mu_k^{(p)} v_k^{(p)}(\mathbf{x}_1) d\mathbf{x}_1 \int_{\partial\Omega} \mu_k^{(p)} [v_k^{(p)}(\mathbf{x}_2)]^* d\mathbf{x}_2 \\ & \times \underbrace{\int_{\Omega} \tilde{G}_0(\mathbf{x}_1, p|\mathbf{x}) \tilde{G}_0(\mathbf{x}_2, p|\mathbf{x}) d\mathbf{x}}_{=-\partial_p \tilde{G}_0(\mathbf{x}_2, p|\mathbf{x}_1)} \\ &= -[\mu_k^{(p)}]^2 \int_{\partial\Omega} v_k^{(p)}(\mathbf{x}_1) d\mathbf{x}_1 \int_{\partial\Omega} [v_k^{(p)}(\mathbf{x}_2)]^* \partial_p \left( \sum_{j=0}^{\infty} \frac{[v_j^{(p)}(\mathbf{x}_1)]^* v_j^{(p)}(\mathbf{x}_2)}{\mu_j^{(p)}} \right) d\mathbf{x}_2, \end{aligned} $$

where we substituted the expansion (31). Writing

$$ \partial_p \left( \frac{[v_j^{(p)}(\mathbf{x}_1)]^* v_j^{(p)}(\mathbf{x}_2)}{\mu_j^{(p)}} \right) = \frac{\partial_p [v_j^{(p)}(\mathbf{x}_1)]^* v_j^{(p)}(\mathbf{x}_2)}{\mu_j^{(p)}} + \frac{[v_j^{(p)}(\mathbf{x}_1)]^* \partial_p v_j^{(p)}(\mathbf{x}_2)}{\mu_j^{(p)}} - \frac{[v_j^{(p)}(\mathbf{x}_1)]^* v_j^{(p)}(\mathbf{x}_2)}{[\mu_j^{(p)}]^2} \partial_p \mu_k^{(p)}, $$

one can separately evaluate three contributions by using the orthonormality of eigenfunctions $\{v_k^{(p)}\}$ on $\partial\Omega$:

$$ \begin{aligned} \|V_k^{(p)}\|_{L^2(\Omega)}^2 &= -[\mu_k^{(p)}]^2 \left\{ \left( \int_{\partial\Omega} v_k^{(p)}(\mathbf{x}_1) \frac{\partial_p [v_k^{(p)}(\mathbf{x}_1)]^*}{\mu_k^{(p)}} d\mathbf{x}_1 \right) \right. \\ & \qquad \left. + \left( \int_{\partial\Omega} [v_k^{(p)}(\mathbf{x}_2)]^* \frac{\partial_p v_k^{(p)}(\mathbf{x}_2)}{\mu_k^{(p)}} d\mathbf{x}_2 \right) - \frac{\partial_p \mu_k^{(p)}}{[\mu_k^{(p)}]^2} \right\} \\ &= -\mu_k^{(p)} \underbrace{\left( \int_{\partial\Omega} \partial_p |v_k^{(p)}(\mathbf{x})|^2 d\mathbf{x} \right)}_{=0} + \partial_p \mu_k^{(p)}. \end{aligned} $$
Exchanging the order of integration over $x$ and differentiation with respect to $p$ yields
$\partial_p \|v_k^{(p)}\|_{L^2(\partial\Omega)}^2 = 0$ and thus implies the identity (58).

References

[1] M. E. Taylor, Partial differential equations: basic theory, Vol. 1, Springer Science & Business Media, 1996.

[2] M. S. Agranovich, On a mixed Poincaré-Steklov type spectral problem in a Lipschitz domain, Russian Journal of Mathematical Physics 13 (3) (2006) 239–244.

[3] A. Girouard, I. Polterovich, Spectral geometry of the Steklov problem (survey article), Journal of Spectral Theory 7 (2) (2017) 321–359.

[4] M. Levitin, D. Mangoubi, I. Polterovich, Topics in Spectral Geometry, to appear in the AMS Graduate Studies in Mathematics series, 2023.

[5] M. Cheney, D. Isaacson, J. C. Newell, Electrical impedance tomography, SIAM review 41 (1) (1999) 85–101.

[6] L. Borcea, Electrical impedance tomography, Inverse problems 18 (6) (2002) R99.

[7] M. S. Zhdanov, G. V. Keller, The geoelectrical methods in geophysical exploration, Methods in geochemistry and geophysics 31 (1994) I–IX.

[8] G. Uhlmann, Inverse problems: seeing the unseen, Bulletin of Mathematical Sciences 4 (2014) 209–279.

[9] D. S. Grebenkov, Spectral theory of imperfect diffusion-controlled reactions on heterogeneous catalytic surfaces, The Journal of Chemical Physics 151 (10) (2019) 104108.

[10] D. S. Grebenkov, Paradigm shift in diffusion-mediated surface phenomena, Physical Review Letters 125 (7) (2020) 078102.

[11] D. S. Grebenkov, Surface hopping propagator: An alternative approach to diffusion-influenced reactions, Physical Review E 102 (3) (2020) 032125.

[12] P. D. Hislop, C. V. Lutzer, Spectral asymptotics of the Dirichlet-to-Neumann map on multiply connected domains in rd, Inverse Problems 17 (6) (2001) 1717.

[13] W. Arendt, A. F. ter Elst, The Dirichlet-to-Neumann operator on rough domains, Journal of Differential Equations 251 (8) (2011) 2100–2124.

[14] W. Arendt, A. F. ter Elst, The Dirichlet-to-Neumann operator on exterior domains, Potential Anal 43 (2) (2015) 313–340.

[15] I. Polterovich, D. A. Sher, J. A. Toth, Nodal length of Steklov eigenfunctions on real-analytic Riemannian surfaces, Journal für die reine und angewandte Mathematik (Crelles Journal) 2019 (754) (2019) 17–47.
[16] J. Galkowski, J. A. Toth, Pointwise bounds for Steklov eigenfunctions, The Journal of Geometric Analysis 29 (2019) 142–193.

[17] T. Daudé, B. Helffer, F. Nicoleau, Exponential localization of Steklov eigenfunctions on warped product manifolds: the flea on the elephant phenomenon, Annales mathématiques du Québec (2021) 1–36.

[18] B. Helffer, A. Kachmar, Semi-classical edge states for the Robin Laplacian, Matematika 68 (2) (2022) 454–485.

[19] M. Levitin, L. Parnovski, I. Polterovich, D. A. Sher, Sloshing, Steklov and corners: Asymptotics of Steklov eigenvalues for curvilinear polygons, Proceedings of the London Mathematical Society 125 (3) (2022) 359–487.

[20] A. Girouard, M. Karpukhin, M. Levitin, I. Polterovich, The Dirichlet-to-Neumann map, the boundary Laplacian, and Hörmander's rediscovered manuscript, Journal of Spectral Theory 12 (1) (2022) 195–225.

[21] B. Colbois, A. Girouard, C. Gordon, D. Sher, Some recent developments on the steklov eigenvalue problem, Revista Matemática Complutense (2023) 1–161.

[22] H. Brezis, Functional analysis, Sobolev spaces and partial differential equations, Vol. 2, Springer, 2011.

[23] W. Stekloff, Sur les problemes fondamentaux de la physique mathématique (suite et fin), in: Annales Scientifiques de l'École Normale Supérieure, Vol. 19, 1902, pp. 455–490.

[24] N. Kuznetsov, T. Kulczycki, M. Kwaśnicki, A. Nazarov, S. Poborchi, I. Polterovich, B. Siudeja, The legacy of Vladimir Andreevich Steklov, Notices of the AMS 61 (1) (2014) 190.

[25] A. B. Andreev, T. D. Todorov, Isoparametric finite-element approximation of a Steklov eigenvalue problem, IMA journal of numerical analysis 24 (2) (2004) 309–322.

[26] H. Bi, Y. Yang, A two-grid method of the non-conforming Crouzeix–Raviart element for the Steklov eigenvalue problem, Applied Mathematics and Computation 217 (23) (2011) 9669–9678.

[27] Q. Li, Y. Yang, A two-grid discretization scheme for the Steklov eigenvalue problem, Journal of Applied Mathematics and Computing 36 (1-2) (2011) 129–139.

[28] H. Bi, H. Li, Y. Yang, An adaptive algorithm based on the shifted inverse iteration for the Steklov eigenvalue problem, Applied Numerical Mathematics 105 (2016) 64–81.

[29] H. Xie, A type of multilevel method for the Steklov eigenvalue problem, IMA Journal of Numerical Analysis 34 (2) (2014) 592–608.
[30] Q. Li, Q. Lin, H. Xie, Nonconforming finite element approximations of the Steklov eigenvalue problem and its lower bound approximations, *Applications of Mathematics* 58 (2) (2013) 129–151.

[31] Y. Yang, Q. Li, S. Li, Nonconforming finite element approximations of the Steklov eigenvalue problem, *Applied Numerical Mathematics* 59 (10) (2009) 2388–2401.

[32] E. Akhmetgaliyev, C.-Y. Kao, B. Osting, Computational methods for extremal Steklov problems, *SIAM Journal on Control and Optimization* 55 (2) (2017) 1226–1240.

[33] O. P. Bruno, J. Galkowski, Domains without dense Steklov nodal sets, *Journal of Fourier Analysis and Applications* 26 (3) (2020) 45.

[34] J.-T. Chen, J.-W. Lee, K.-T. Lien, Analytical and numerical studies for solving Steklov eigenproblems by using the boundary integral equation method/boundary element method, *Engineering Analysis with Boundary Elements* 114 (2020) 136–147.

[35] V. D. Kupradze, M. A. Aleksidze, The method of functional equations for the approximate solution of certain boundary value problems, *USSR Computational Mathematics and Mathematical Physics* 4 (4) (1964) 82–126.

[36] B. Bogosel, The method of fundamental solutions applied to boundary eigenvalue problems, *Journal of Computational and Applied Mathematics* 306 (2016) 265–285.

[37] W. Alhejaili, C.-Y. Kao, Numerical studies of the Steklov eigenvalue problem via conformal mappings, *Applied Mathematics and Computation* 347 (2019) 785–802.

[38] D. S. Grebenkov, Encounter-based approach to the escape problem, *Physical Review E* 107 (4) (2023) 044105.

[39] A. A. Lacey, J. Sabina, J. R. Ockendon, Multidimensional reaction diffusion equations with nonlinear boundary conditions, *SIAM Journal on Applied Mathematics* 58 (5) (1998) 1622–1647.

[40] M. Levitin, L. Parnovski, On the principal eigenvalue of a Robin problem with a large parameter, *Mathematische Nachrichten* 281 (2) (2008) 272–281.

[41] M. Khalile, Spectral asymptotics for Robin Laplacians on polygonal domains, *Journal of Mathematical Analysis and Applications* 461 (2) (2018) 1498–1543.

[42] M. Khalile, K. Pankrashkin, Eigenvalues of Robin Laplacians in infinite sectors, *Mathematische Nachrichten* 291 (5-6) (2018) 928–965.

[43] R. S. Laugesen, The Robin Laplacian—spectral conjectures, rectangular theorems, *Journal of Mathematical Physics* 60 (12) (2019).