# The hot spots conjecture can be false: Some numerical examples

Andreas Kleefeld1111Author to whom any correspondence should be addressed. 1
Forschungszentrum Jülich GmbH, Jülich Supercomputing Centre, 52425 Jülich,
Germany<EMAIL_ADDRESS>
###### Abstract

The hot spots conjecture is only known to be true for special geometries. It
can be shown numerically that the hot spots conjecture can fail to be true for
easy to construct bounded domains with one hole. The underlying eigenvalue
problem for the Laplace equation with Neumann boundary condition is solved
with boundary integral equations yielding a non-linear eigenvalue problem. Its
discretization via the boundary element collocation method in combination with
the algorithm by Beyn yields highly accurate results both for the first non-
zero eigenvalue and its corresponding eigenfunction which is due to
superconvergence. Additionally, it can be shown numerically that the ratio
between the maximal/minimal value inside the domain and its maximal/minimal
value on the boundary can be larger than $1+10^{-3}$. Finally, numerical
examples for easy to construct domains with up to five holes are provided
which fail the hot spots conjecture as well.

###### ams:

35J25, 35P20, 65F15, 65M38, 78A46

††: arXiv

  * 4 January 2021

Keywords: interior Neumann eigenvalues, Helmholtz equation, potential theory,
boundary integral equations, numerics

##  1 Introduction

The hot spots conjecture has been given in 1974 by Jeffrey Rauch [32] and
explicitly stated a decade later in Kawohl [19]. Refer also to the paper by
Bañuelos & Burdzy [6] from 1999. Since then a lot of researchers have worked
on this challenging problem (see Judge & Mondal [18] for a recent overview
from 2020).

Before stating it in mathematical terms, we explain it in a simple fashion.
Imagine that we have a flat piece of metal $D$ where $D$ is a bounded subset
of the two-dimensional space (Euclidean domain) which can have holes with a
sufficiently smooth boundary. Next, an (almost) arbitrary initial temperature
distribution is provided on $D$ (refer to [6, p. 2]). Assume that the domain
is insulated, then the hottest and coldest spot of $D$ will appear on the
boundary when waiting for a long time.

Now, we go into the mathematical detail: That means we have to solve the heat
equation $\partial_{t}u=\Delta u$ for $t\rightarrow\infty$ with homogeneous
Neumann boundary condition $\partial_{\nu}u=0$ and ‘almost’ arbitrary initial
condition in an open connected bounded $D$ with Lipschitz boundary for its
equilibrium (see [6, p. 2] and [34] for the definition of a Lipschitz domain).
Refer to Figure 1 for an example.

(a) Solution at time $t_{1}=1/200$ (b) Solution at time $t_{2}=1/10$ (c)
Solution at time $t_{3}=1/2$ (d) Solution at time $t_{4}=2$

Figure 1:  The numerical solution of the heat equation
$\partial_{t}u=\frac{1}{10}\Delta u$ and initial condition with standard
normal random numbers and homogeneous Neumann boundary condition for time
$t_{1}=1/200$, $t_{2}=1/10$, $t_{3}=1/2$, and $t_{4}=2$ using $h=1/100$ and
$k=1/100$ (see also [2, p. 11] for more details on the implementation of the
exponential time differencing method). The maximal and minimal value for
$T_{4}$ appear on the boundary. Note that the solution $T_{4}$ is
approximately representing the first non-zero Neumann eigenfunction (see also
Figure 9).

Precisely, we have to find the smallest non-trivial eigenvalue of the
Laplacian with homogeneous Neumann boundary condition and the corresponding
eigenfunction. Note that the smallest eigenvalue of the Laplacian with
homogeneous Neumann boundary condition is zero with corresponding
eigenfunction $u_{0}=\mathrm{const}$. Mathematically, we have to find a
solution $u\neq 0$ and the smallest $k\in\mathbb{R}_{>0}$ such that the
Helmholtz equation $\Delta u+k^{2}u=0$ in $D$ with $\partial_{\nu}u=0$ on the
boundary $\Gamma$. All solutions $k$ are called non-trivial interior Neumann
eigenvalues and $\lambda_{i}=k^{2}_{i}$ will be the $i$-th non-trivial Neumann
eigenvalue of the Laplacian. Its corresponding eigenfunctions are denoted by
$u_{i}$. Further, it is known that the eigenvalues satisfy
$0=\lambda_{0}<\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}\leq\ldots$ when $D$
is a bounded planar domain with Lipschitz boundary (see for example [16, p.
449] and the references therein, specifically [33]). If
$\lambda_{1}<\lambda_{2}$ and $\left\langle u(0,x),u_{1}\right\rangle\neq 0$,
then $u(t,x)=\mathrm{e}^{\lambda_{1}t}\left\langle u(0,x),u_{1}\right\rangle+$
lower terms. Note that the first non-trivial eigenvalue can have multiplicity
more than one which means that there can be more than one eigenfunction.

Now, the conjecture can be stated as (refer also to [6, p. 2]): Let
$D\subset\mathbb{R}^{2}$ be an open connected bounded domain with Lipschitz
boundary $\Gamma$. Then:

  * C1:

For each eigenfunction $u_{2}(x)$ corresponding to $\lambda_{2}$ which is not
identically zero, we have

$\inf_{x\in\Gamma}u_{2}(x)<u_{2}(y)<\sup_{x\in\Gamma}u_{2}(x)\quad\forall y\in
D\,.$

  * C2:

For each eigenfunction $u_{2}(x)$ corresponding to $\lambda_{2}$ which is not
identically zero, we have

$\inf_{x\in\Gamma}u_{2}(x)\leq
u_{2}(y)\leq\sup_{x\in\Gamma}u_{2}(x)\quad\forall y\in D\,.$

  * C3:

There exist an eigenfunction $u_{2}(x)$ corresponding to $\lambda_{2}$ which
is not identically zero, such that

$\inf_{x\in\Gamma}u_{2}(x)\leq
u_{2}(y)\leq\sup_{x\in\Gamma}u_{2}(x)\quad\forall y\in D\,.$

Here, C1 is the original conjecture of Rauch. The hypothesis has been shown to
be true for some special geometries such as parallelepipeds, balls,
rectangles, cylinders [19], obtuse triangles [6], some convex and non-convex
domains with symmetry [6], wedges [3], lip domains [4], convex domains with
two axes of symmetry [17], convex $C^{1,\alpha}$ domains ($0<\alpha<1$) with
one axis of symmetry [31], a certain class of planar convex domains [29],
subequilateral isosceles triangles [30], a certain class of acute triangles
[37], Euclidean triangles [18], and strips on two-dimensional Riemannian
manifolds [25].

It is assumed that the hot spots conjecture is true for arbitrary convex
domains, but a proof is still open. The hot spots conjecture is assumed to be
true also for simply-connected bounded non-convex domains, but no successful
attempts (neither theoretically nor numerically) have been made to prove this
conjecture or to find a counterexample.

It has been shown that for some domains with one hole that the hot spots
conjecture is true (for example an annulus [19]), but that there are also
domains with one or more holes where the hot spots conjecture is false (see
Burdzy [9], Burdzy & Werner [10], and Bass & Burdzy [7], respectively). For
domains on manifolds, we refer the reader to [13].

However, the proofs in [9, 10, 7] are very technical and are based on
stochastic arguments. No numerical results support their counterexamples,
since their domains are too complicated for being constructed. To be precise,
they are very thin and have a polygonal structure. Further, the first non-
trivial Neumann eigenvalue is assumed to be simple. If this is not the case,
the proof collapses.

The only non-published numerical results given so far are for triangles in the
PolyMath project 7 ‘Hot spots conjecture’ from 2012 to 2013 using the finite
element method.

Further work in the direction of better understanding the conjecture is given
for example by Steinerberger [39]. Related results on graphs are given by
Lederman & Steinerberger [27].

## Contribution

It is the goal of this paper to construct ‘simple’ domains with one hole and
show numerically with high precision (due to superconvergence) that those
domains do not satisfy the hot spots conjecture. The method based on boundary
integral equations is very efficient and its convergence is faster than
expected. Additionally, we show the influence on the location of the hot spots
by changing the boundary of the domain in order to understand this connection.
It is believed that this might help researchers to provide assumptions when
the hot spots conjecture will be true or false for arbitrary bounded simply-
connected domains which are not necessarily convex. We show that it is
possible to construct domains with one hole such that the ratio between the
maximum/minimum in the interior and its maximum/minimum on the boundary is
larger than $1+10^{-3}$. Finally, numerical results are given that show that
there exist domains with up to five holes which do not satisfy the hot spots
conjecture as well. The Matlab programs including the produced data are
available at github https://github.com/kleefeld80/hotspots and can be used by
any researcher trying their own geometries and to reproduce the numerical
results within this article.

## Outline of the paper

In Section 2, we explain the algorithm in order to compute the first non-zero
Neumann Laplace eigenvalue and its corresponding eigenfunction for an
arbitrary domain with or without a hole using boundary integral equations
resulting in a non-linear eigenvalue problem. Further, it is shown in detail
how to discretize the boundary integral equations via the boundary element
collocation method and how to numerically solve the non-linear eigenvalue
problem. Extensive numerical results are provided in Section 3 showing the
superconvergence and highly accurate results for domains with one or no hole.
Domains are provided that show the failure of the hot spots conjecture and
further interesting results. The extension to domains with up to five holes is
straightforward and given at the end of this section as well. A short summary
and outlook is given in Section 4.

##  2 The algorithm

In this section, we explain the algorithm to compute numerically non-trivial
interior Neumann eigenvalues and its corresponding eigenfunction to high
accuracy for bounded domains with one hole (the extension to more than one
hole is straightforward) very efficiently. The ingredients are boundary
integral equations and its approximation via boundary element collocation
method; that is, a two-dimensional problem is reduced to a one-dimensional
problem. The resulting non-linear eigenvalue problem is solved using complex-
valued contour integrals integrating over the resolvent reducing the non-
linear eigenvalue problem to a linear eigenvalue problem which is possible due
to Keldysh’s theorem (see Beyn [8]).

###  2.1 Notations

We consider a bounded Lipschitz domain $D\subset\mathbb{R}^{2}$ with one hole.
The outer boundary $\Gamma_{1}$ is assumed to be sufficiently smooth that is
oriented counter-clockwise and a sufficiently smooth inner boundary
$\Gamma_{2}$ that is oriented clockwise. The normal $\nu_{1}$ on the boundary
$\Gamma_{1}$ is pointing into the unbounded exterior $E$. The normal $\nu_{2}$
on the boundary $\Gamma_{2}$ is pointing into the bounded exterior $I$. We
refer the reader to Figure 2. The boundary of $D$ is given by
$\Gamma=\Gamma_{1}\cup\Gamma_{2}$ ($\Gamma_{1}\cap\Gamma_{2}\neq\emptyset$).

$D$$E$$\Gamma_{1}$$\nu_{1}$$I$$\Gamma_{2}$$\nu_{2}$ Figure 2:  Used notations
for a bounded domain with one hole.

Note that we also consider bounded domains without a hole. In this case, we
have $\Gamma_{2}=\emptyset$ and hence $\Gamma=\Gamma_{1}$ and $I=\emptyset$.

###  2.2 Boundary integral equation

The solution to the Helmholtz equation $\Delta u+k^{2}u=0$ (reduced wave
equation) in the domain $D$ for a given wave number $k$ with
$\mathrm{Im}(k)\geq 0$ is given by (see [12, Theorem 2.1])

$\displaystyle
u(x)=\int_{\Gamma}\partial_{\nu(y)}u(y)\cdotp\Phi_{k}(x,y)-u(y)\cdotp\partial_{\nu(y)}\Phi_{k}(x,y)\,\mathrm{d}s(y)\,,\quad
x\in D$

which can be written in our notation as

$\displaystyle u(x)$ $\displaystyle=$
$\displaystyle\int_{\Gamma_{1}}\partial_{\nu_{1}(y)}u(y)\cdotp\Phi_{k}(x,y)-u(y)\cdotp\partial_{\nu_{1}(y)}\Phi_{k}(x,y)\,\mathrm{d}s(y)$
(1) $\displaystyle+$
$\displaystyle\int_{\Gamma_{2}}\partial_{\nu_{2}(y)}u(y)\cdotp\Phi_{k}(x,y)-u(y)\cdotp\partial_{\nu_{2}(y)}\Phi_{k}(x,y)\,\mathrm{d}s(y)\,,\quad
x\in D$

where $\Phi_{k}(x,y)=\mathrm{i}H_{0}^{(1)}(k\|x-y\|)/4$, $x\neq y$ denotes the
fundamental solution of the Helmholtz equation in two dimensions (see [12, p.
66]). Here, $H_{0}^{(1)}$ denotes the first-kind Hankel function of order
zero. We denote $u(y)$ for $y\in\Gamma_{1}$ as $u_{1}(y)$ and similarly $u(y)$
for $y\in\Gamma_{2}$ as $u_{2}(y)$. Hence, we can write (1) as

$\displaystyle u(x)$ $\displaystyle=$
$\displaystyle-\int_{\Gamma_{1}}u_{1}(y)\cdotp\partial_{\nu_{1}(y)}\Phi_{k}(x,y)\,\mathrm{d}s(y)$
(2)
$\displaystyle-\int_{\Gamma_{2}}u_{2}(y)\cdotp\partial_{\nu_{2}(y)}\Phi_{k}(x,y)\,\mathrm{d}s(y)\,,\quad
x\in D$

where we also used the homogeneous Neumann boundary conditions
$\partial_{\nu_{1}}u=0$ and $\partial_{\nu_{2}}u=0$. We rewrite (2) as

$\displaystyle
u(x)=-\mathrm{DL}_{k}^{\Gamma_{1}}u_{1}(x)-\mathrm{DL}_{k}^{\Gamma_{2}}u_{2}(x)\,,\quad
x\in D$ (3)

where we used the notation

$\mathrm{DL}_{k}^{\Gamma_{i}}\psi_{i}(x)=\int_{\Gamma_{i}}\psi_{i}(y)\cdotp\partial_{\nu_{i}(y)}\Phi_{k}(x,y)\,\mathrm{d}s(y)\,,\quad
x\in D\,,\quad i=1,2$

for the acoustic double layer potential with density $\psi_{i}$ (see [12, p.
39]). Assume for a moment that $k$ is given. The functions $u_{1}$ and $u_{2}$
are still unknown. Once we know them, we can compute the solution $u$ inside
the domain of $D$ at any point we want using (3).

Now, we explain how to obtain those functions $u_{1}$ and $u_{2}$ on the
boundary. Letting $x\in D$ approach the boundary $\Gamma_{1}$ and using the
jump relation of the acoustic double layer operator (see [12, p. 39] for the
smooth boundary case, otherwise [34]), yields the boundary integral equation

$\displaystyle
u_{1}(x)=-\left(\mathrm{D}_{k}^{\Gamma_{1}\rightarrow\Gamma_{1}}u_{1}(x)-\left(1-\Omega_{1}(x)\right)u_{1}(x)\right)-\mathrm{D}_{k}^{\Gamma_{2}\rightarrow\Gamma_{1}}u_{2}(x)\,,\quad
x\in\Gamma_{1}$ (4)

where we used the notation

$\mathrm{D}_{k}^{\Gamma_{i}\rightarrow\Gamma_{j}}\psi_{i}(x)=\int_{\Gamma_{i}}\psi_{i}(y)\cdotp\partial_{\nu_{i}(y)}\Phi_{k}(x,y)\,\mathrm{d}s(y)\,,\quad
x\in\Gamma_{j}\,,\quad i,j=1,2$

for the double layer operator (see [12, p. 41]). Here, $\Omega_{1}(x)$ denotes
the interior solid angle at a point $x$ on $\Gamma_{1}$. When the boundary is
smooth at this point, then $\Omega_{1}(x)=1/2$. In fact, $\Omega_{1}(x)$ is
$1/2$ almost everywhere for Lipschitz domains.

Similarly, we obtain for $x\in D$ approaching the boundary $\Gamma_{2}$ and
using the jump relation for the double layer operator

$\displaystyle
u_{2}(x)=-\mathrm{D}_{k}^{\Gamma_{1}\rightarrow\Gamma_{2}}u_{1}(x)-\left(\mathrm{D}_{k}^{\Gamma_{2}\rightarrow\Gamma_{2}}u_{2}(x)-\left(1-\Omega_{2}(x)\right)u_{2}(x)\right)\,,\quad
x\in\Gamma_{2}\,.$ (5)

We can rewrite (4) and (5) as a $2\times 2$ system of boundary integral
equations in the form

$\displaystyle\underbrace{\left(\overbrace{\left(\begin{matrix}\Omega_{1}&0\\\
0&\Omega_{2}\end{matrix}\right)}^{\mathrm{C}}\overbrace{\left(\begin{matrix}\mathrm{I}&0\\\
0&\mathrm{I}\end{matrix}\right)}^{\mathrm{I}_{\mathrm{B}}}+\overbrace{\left(\begin{matrix}\mathrm{D}_{k}^{\Gamma_{1}\rightarrow\Gamma_{1}}&\mathrm{D}_{k}^{\Gamma_{2}\rightarrow\Gamma_{1}}\\\
\mathrm{D}_{k}^{\Gamma_{1}\rightarrow\Gamma_{2}}&\mathrm{D}_{k}^{\Gamma_{2}\rightarrow\Gamma_{2}}\end{matrix}\right)}^{\mathrm{K}(k)}\right)}_{\mathrm{M}(k)}\underbrace{\left(\begin{matrix}u_{1}\\\
u_{2}\end{matrix}\right)}_{u}=\left(\begin{matrix}0\\\
0\end{matrix}\right)\quad\text{on }\Gamma$ (6)

where $\mathrm{I}$ and $\mathrm{I}_{\mathrm{B}}$ denotes the identity and the
$2\times 2$ block identity operator, respectively. Hence, we have to
numerically solve the non-linear eigenvalue problem (6) written as
$\mathrm{M}(k)u=0$ to find the smallest non-trivial (real) eigenvalue $k$ and
the corresponding eigenfunction $u$. Then, we can numerically evaluate (3) to
compute the eigenfunction at any point in the interior we want. As in [21, p.
188], we can argue that the compact operator $\mathrm{K}(k)$ maps from
$\mathcal{H}^{-1/2}(\Gamma_{1})\times\mathcal{H}^{-1/2}(\Gamma_{2})$ to
$\mathcal{H}^{1/2}(\Gamma_{1})\times\mathcal{H}^{1/2}(\Gamma_{2})$. Here,
$\mathcal{H}^{s}(\Gamma)$ denotes a Sobolev space of order $s\in\mathbb{R}$ on
the domain $\Gamma$ which are defined via Bessel potentials (see [28, pp.
75–76] for more details). The operator
$\mathrm{M}(k)=\mathrm{C}\cdotp\mathrm{I}_{\mathrm{B}}+\mathrm{K}(k)$ is
Fredholm of index zero for $k\in\mathbb{C}\backslash\mathbb{R}_{\leq 0}$ and
therefore the theory of eigenvalue problems for holomorphic Fredholm operator-
valued functions applies to $\mathrm{M}(k)$.

###  2.3 Discretization

In this section, we explain how to discretize (6) using quadratic
interpolation of the boundary, but using piecewise quadratic interpolation
with $\alpha=(1-\sqrt{3/5})/2$ (see [22] for the 3D case) instead of quadratic
interpolation for the unknown $u$ on each of the $n_{f}$ boundary elements
which ultimately leads to the non-linear eigenvalue

$\displaystyle\mathbf{M}(k)\vec{u}=\vec{0}$ (7)

where the matrix is of size $3\cdotp 2\cdotp n_{f}\times 3\cdotp 2\cdotp
n_{f}$. The size of the matrix is slightly larger than the one given in
Kleefeld [21], but it has the advantage that no singular integral has to be
evaluated numerically, since we can use a similar singularity subtraction
technique as explained in Kleefeld and Lin [23, pp. A1720–A1721] and the
convergence rate is slightly higher. The details are about to follow for a
domain without a hole for simplicity. In this case, we have to solve a
boundary integral equation of the second kind of the form

$\displaystyle\Omega_{1}(x)u_{1}(x)+\int_{\Gamma_{1}}u_{1}(y)\partial_{\nu_{1}(y)}\Phi_{k}(x,y)\,\mathrm{d}s(y)=0\,,\quad
x\in\Gamma_{1}\,.$ (8)

First, we subdivide the boundary $\Gamma_{1}$ into $n_{f}$ pieces denoted by
$\Delta_{j}$ with $j=1,\ldots,n_{f}$. A subdivision into four pieces is shown
in Figure 3 for the unit circle.

$\Delta_{1}$$\Delta_{2}$$\Delta_{3}$$\Delta_{4}$ Figure 3: Subdivision of the
given boundary $\Gamma_{1}$ into four pieces $\Delta_{1}$, $\Delta_{2}$,
$\Delta_{3}$, and $\Delta_{4}$.

Then equation (8) can be equivalently written as

$\displaystyle\Omega_{1}(x)u_{1}(x)+\sum_{j=1}^{n_{f}}\int_{\Delta_{j}}u_{1}(y)\partial_{\nu_{1}(y)}\Phi_{k}(x,y)\,\mathrm{d}s(y)=0\,,\quad
x\in\Gamma_{1}\,.$

For each $j$ there exists a unique map $m_{j}$ which maps from the standard
interval $\sigma=[0,1]$ to $\Delta_{j}$. Then, we can apply a simple change of
variables to each integral over $\Delta_{j}$ giving

$\displaystyle\Omega_{1}(x)u_{1}(x)$ $\displaystyle+$
$\displaystyle\sum_{j=1}^{n_{f}}\int_{\sigma}u_{1}(m_{j}(s))\partial_{\nu_{1}(m_{j}(s))}\Phi_{k}(x,m_{j}(s))J(s)\,\mathrm{d}s(s)=0\,,\;x\in\Gamma_{1}$

with the Jacobian given by $J(s)=\|\partial_{s}m_{j}(s)\|$. In most cases, we
can explicitly write down this map. However, we approximate each $m_{j}(s)$ by
a quadratic interpolation polynomial
$m_{j}(s)\approx\widetilde{m}_{j}(s)=\sum_{i=1}^{3}v_{(2i-j)\,\mathrm{mod}\,(2n_{f})}L_{i}(s)$
where the Lagrange basis functions are

$L_{1}(s)=u\cdotp(1-2s)\,,\quad L_{2}(s)=4s\cdotp u\,,\quad\text{ and }\quad
L_{3}(s)=s\cdotp(2s-1)$

with $u=1-s$. Here,
$\gamma\,\mathrm{mod}\,\delta=\gamma-\lfloor\frac{\gamma}{\delta}\rfloor\cdotp\delta$
with $\lfloor\cdotp\rfloor$ the floor function. The nodes $v_{\ell}$
($\ell=1,\ldots,2n_{f}$) are the given vertices and midpoints of the $n_{f}$
faces. We refer the reader to Figure 4 for an example with four faces and
eight nodes (four vertices and midpoints, respectively).

$v_{2}$$v_{8}$$v_{6}$$v_{4}$$v_{1}$$v_{7}$$v_{5}$$v_{3}$ Figure 4: The eight
nodes $v_{1},\ldots,v_{8}$, the four vertices $v_{1}$, $v_{3}$, $v_{5}$, and
$v_{7}$ (marked with $\times$), and the four midpoints $v_{2}$, $v_{4}$,
$v_{6}$, and $v_{8}$ (marked with $\circ$) for the four faces $\Delta_{1}$,
$\Delta_{2}$, $\Delta_{3}$, and $\Delta_{4}$.

An example how the approximation of $\Delta_{1}$ via a quadratic interpolation
polynomial using two vertices and the midpoint looks like is shown in Figure
5.

$\frac{\pi}{4}$$0$$\frac{\pi}{2}$ Figure 5: The approximation of the first
part of the boundary $\Delta_{1}$ (solid line) via a quadratic interpolation
polynomial using two vertices and the midpoint (dashed line).

Next, we define the ‘collocation nodes’ $\widetilde{v}_{j,k}$ by
$\widetilde{v}_{j,k}=\widetilde{m}_{j}(q_{k})$ for $j=1,\ldots,n_{f}$ and for
$k=1,2,3$ where $q_{1}=\alpha$, $q_{2}=1/2$, and $q_{3}=1-\alpha$ with
$0<\alpha<1/2$ a given and fixed constant. This ensures that the collocation
nodes are always lying within a piece of the boundary and at those points the
interior solid angle is $1/2$. For a specific choice of $\alpha$ the overall
convergence rate can be improved. The first three collocation nodes on the
approximated boundary for the unit circle using $\alpha=(1-\sqrt{3/5})/2$ are
shown in Figure 6.

Figure 6: The first three collocation nodes (solid triangles) on the first
part of the approximated boundary (dashed line) for the unit circle using
$\alpha=(1-\sqrt{3/5})/2$. The exact boundary is shown with a solid line
including the two vertices marked by a cross.

The unknown function $u_{1}(\widetilde{m}_{j}(s))$ is now approximated on each
of the $j$ pieces by a quadratic interpolation polynomial of the form
$\sum_{k=1}^{3}u_{1}(\widetilde{m}_{j}(q_{k}))\widetilde{L}_{k}(s)$ which can
be written as $\sum_{k=1}^{3}u_{1}(\widetilde{v}_{j,k})\widetilde{L}_{k}(s)$
where the Lagrange basis functions are given by

$\displaystyle\widetilde{L}_{1}(s)=\frac{u-\alpha}{1-2\alpha}\frac{1-2s}{1-2\alpha}\,,\;\widetilde{L}_{2}(s)=4\frac{s-\alpha}{1-2\alpha}\frac{u-\alpha}{1-2\alpha}\,,\;\widetilde{L}_{3}(s)=\frac{s-\alpha}{1-2\alpha}\frac{2s-1}{1-2\alpha}$

with $u=1-s$. We obtain

$\displaystyle\Omega_{1}(x)u_{1}(x)$ $\displaystyle+$
$\displaystyle\sum_{j=1}^{n_{f}}\sum_{k=1}^{3}\int_{\sigma}\partial_{\nu_{1}(\widetilde{m}_{j}(s))}\Phi_{k}(x,\widetilde{m}_{j}(s))\|\partial_{s}\widetilde{m}_{j}(s)\|\widetilde{L}_{k}(s)\,\mathrm{d}s(s)u_{1}(\widetilde{v}_{j,k})=r(x)$

with $r(x)$ the residue which is due to the different approximations. We set
$r(\widetilde{v}_{i,\ell})=0$ and since
$\Omega_{1}(\widetilde{v}_{i,\ell})=1/2$ always by the choice of the
collocation nodes, we obtain the linear system of size $3n_{f}\times 3n_{f}$

$\displaystyle\frac{1}{2}u_{1}(\widetilde{v}_{i,\ell})+\sum_{j=1}^{n_{f}}\sum_{k=1}^{3}a_{i,\ell,j,k}u_{1}(\widetilde{v}_{j,k})=0$

with the resulting integrals

$\displaystyle
a_{i,\ell,j,k}=\int_{\sigma}\partial_{\nu_{1}(\widetilde{m}_{j}(s))}\Phi_{k}(\widetilde{v}_{i,\ell},\widetilde{m}_{j}(s))\|\partial_{s}\widetilde{m}_{j}(s)\|\widetilde{L}_{k}(s)\,\mathrm{d}s(s)$
(9)

which will be approximated by the adaptive Gauss-Kronrod quadrature (see
[36]). This can be written abstractly as $\mathbf{M}(k)\vec{u}=\vec{0}$. Note
that the integrand of the integral of the right-hand side (9) can easily by
written down as

$\frac{\mathrm{i}kH_{1}^{(1)}(kr)}{4r}\left(a\cdotp n_{1}+b\cdotp
n_{2}\right)\widetilde{L}_{k}(s)$

with $H_{1}^{(1)}$ the first-kind Hankel function of order one and with

$\displaystyle a$ $\displaystyle=$
$\displaystyle\left[\widetilde{v}_{i,\ell}-\widetilde{m}_{j}(s)\right]_{1}\,,\quad
b=\left[\widetilde{v}_{i,\ell}-\widetilde{m}_{j}(s)\right]_{2}\,,\quad
r=\sqrt{a^{2}+b^{2}}\,,$ $\displaystyle n_{1}$ $\displaystyle=$
$\displaystyle\left[\partial_{s}\widetilde{m}_{j}(s)\right]_{2}\,,\quad
n_{2}=-\left[\partial_{s}\widetilde{m}_{j}(s)\right]_{1}\,.$

Note that the Jacobian cancels out. A word has to be spent on the following
issue: When $P\neq Q$ with $P=\widetilde{v}_{i,\ell}$ and
$Q=\widetilde{m}_{j}(s)$, then the integrand of the integral of the right-hand
side (9) is smooth. However, for the case $P=Q$ a singularity within the
integral of the right-hand side (9) is present. In this case, we can use the
singularity subtraction method to rewrite the singular integral in the
following form

$\displaystyle\int_{\sigma}\partial_{\nu_{1}(\widetilde{m}_{j}(s))}\Phi_{k}(\widetilde{v}_{i,\ell},\widetilde{m}_{j}(s))\|\partial_{s}\widetilde{m}_{j}(s)\|\widetilde{L}_{k}(s)\,\mathrm{d}s(s)$
$\displaystyle=$
$\displaystyle\int_{\sigma}\partial_{\nu_{1}(\widetilde{m}_{j}(s))}\left(\Phi_{k}(\widetilde{v}_{i,\ell},\widetilde{m}_{j}(s))-\Phi_{0}(\widetilde{v}_{i,\ell},\widetilde{m}_{j}(s))\right)\|\partial_{s}\widetilde{m}_{j}(s)\|\widetilde{L}_{k}(s)\,\mathrm{d}s(s)$
$\displaystyle+$
$\displaystyle\int_{\sigma}\partial_{\nu_{1}(\widetilde{m}_{j}(s))}\Phi_{0}(\widetilde{v}_{i,\ell},\widetilde{m}_{j}(s))\|\partial_{s}\widetilde{m}_{j}(s)\|\widetilde{L}_{k}(s)\,\mathrm{d}s(s)=I_{i,\ell}^{\text{smooth}}+I_{i,\ell}^{\text{sing}}$

where $\Phi_{0}(P,Q)=-\log(|P-Q|)/(2\pi)$ is the fundamental solution of the
Laplace equation. The integral $I_{i,\ell}^{\text{smooth}}$ has a smooth
kernel (no singularity present) and is converging rapidly to zero when
increasing the number of faces (independent of the wave number $k$). Hence, we
directly set $I_{i,\ell}^{\text{smooth}}=0$. The integral
$I_{i,\ell}^{\text{sing}}$ (a singularity is present) can be rewritten as a
sum of integrals without any singularity. This is due to the fact that we have
$\mathrm{D}_{0}^{\Gamma_{1}\rightarrow\Gamma_{1}}\,1\,(x)=-\Omega_{1}(x)$ mit
$\psi=1$, $\forall x\in\Gamma_{1}$ (see [35, p. 363]) and hence, we
approximately use for all $i,\ell$:

$\sum_{j=1}^{n_{f}}\sum_{k=1}^{3}\int_{\sigma}\partial_{\nu_{1}(\widetilde{m}_{j}(s))}\Phi_{0}(\widetilde{v}_{i,\ell},\widetilde{m}_{j}(s))\|\partial_{s}\widetilde{m}_{j}(s)\|\widetilde{L}_{k}(s)\,\mathrm{d}s(s)\approx-\Omega_{1}(\widetilde{v}_{i,\ell})=-\frac{1}{2}$

(that is, the row sum of the matrix $\mathbf{M}(0)$ obtained from the
discretization of the double layer for the Laplace equation shall be $-1/2$)
and hence, we can compute $I_{i,\ell}^{\text{sing}}$ as

$\displaystyle I_{i,\ell}^{\text{sing}}$ $\displaystyle\approx$
$\displaystyle-\underbrace{\Omega_{1}(\widetilde{v}_{i,\ell})}_{\frac{1}{2}}$
(10) $\displaystyle-\mathop{\sum_{j=1}^{n_{f}}}_{j\neq
i}\mathop{\sum_{k=1}^{3}}_{k\neq\ell}\int_{\sigma}\partial_{\nu_{1}(\widetilde{m}_{j}(s))}\Phi_{0}(\widetilde{v}_{i,\ell},\widetilde{m}_{j}(s))\|\partial_{s}\widetilde{m}_{j}(s)\|\widetilde{L}_{k}(s)\,\mathrm{d}s(s)\,.$

Each of the integrands within the integrals of the right-hand side (10) are
smooth and can be computed with the previously mentioned Gauss-Kronrod
quadrature. Note that we never have to compute the interior solid angle since
it is always $1/2$ by the choice of the collocation points. In fact, we never
have to use the value $1/2$ since it cancels out with the $1/2$ within the
definition of $\mathbf{M}(k)$.

Finally, if we use a domain with a hole, we obtain the system of size $m\times
m=3\cdotp 2\cdotp n_{f}\times 3\cdotp 2\cdotp n_{f}$ when including the
boundary $\Gamma_{2}$ and the unknown function $u_{2}$. Written abstractly we
obtain the non-linear eigenvalue problem

$\displaystyle\mathbf{M}(k)\vec{u}=\vec{0}$ (11)

Of course, the extension to more than one hole is obvious. In this case, the
matrix $\mathbf{M}(k)$ within (11) will be of size $(q-1)\cdotp 2\cdotp
n_{f}\times(q-1)\cdotp 2\cdotp n_{f}$ where $q$ denotes the number of holes
within the domain.

###  2.4 Non-linear eigenvalue problem

The non-linear eigenvalue problem (11) is solved with the Beyn algorithm [8].
It is based on complex-valued contour integrals integrating over the resolvent
reducing the non-linear eigenvalue to a linear one (of very small size) which
can be achieved by Keldysh’s theorem. Precisely, a user-specified
$2\pi$-periodic contour $\gamma$ of class $C^{1}$ within the complex plane has
to be given. We need a contour that is enclosing a part of the real line where
the smallest non-zero eigenvalue is expected. We usually use a circle with
radius $R$ and center $(\mu,0)$ (in order to exclude the eigenvalue zero, we
choose $\mu>R$). In this case, we have
$\varphi(t)=\mu+R\cos(t)+\mathrm{i}R\sin(t)$ which satisfies $\varphi\in
C^{\infty}$. The number of eigenvalues including their multiplicity within the
contour $\gamma$ is denoted by $n(\gamma)$. With the randomly chosen matrix
$\hat{\mathbf{V}}\in\mathbb{C}^{m\times\ell}$ with $m\gg\ell\geq n(\gamma)$
the two contour integrals of the form

$\displaystyle\mathbf{A}_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi\mathrm{i}}\int_{\gamma}\mathbf{M}^{-1}(k)\hat{\mathbf{V}}\,\mathrm{d}s(k)\,,$
$\displaystyle\mathbf{A}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi\mathrm{i}}\int_{\gamma}k\mathbf{M}^{-1}(k)\hat{\mathbf{V}}\,\mathrm{d}s(k)$

over the given contour $\gamma$ are now rewritten as

$\displaystyle\mathbf{A}_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi\mathrm{i}}\int_{0}^{2\pi}\mathbf{M}^{-1}(\varphi(t))\hat{\mathbf{V}}\varphi^{\prime}(t)\,\mathrm{d}s(t)\,,$
$\displaystyle\mathbf{A}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi\mathrm{i}}\int_{0}^{2\pi}\varphi(t)\mathbf{M}^{-1}(\varphi(t))\hat{\mathbf{V}}\varphi^{\prime}(t)\,\mathrm{d}s(t)$

and approximated by the trapezoidal rule yielding

$\displaystyle\mathbf{A}_{0,N}$ $\displaystyle=$
$\displaystyle\frac{1}{\mathrm{i}N}\sum_{j=0}^{N-1}\mathbf{M}^{-1}(\varphi(t_{j}))\hat{\mathbf{V}}\varphi^{\prime}(t_{j}),$
$\displaystyle\mathbf{A}_{1,N}$ $\displaystyle=$
$\displaystyle\frac{1}{\mathrm{i}N}\sum_{j=0}^{N-1}\varphi(t_{j})\mathbf{M}^{-1}(\varphi(t_{j}))\hat{\mathbf{V}}\varphi^{\prime}(t_{j})\,,$

where the parameter $N$ is given and the equidistant nodes are $t_{j}=2\pi
j/N$, $j=0,\ldots,N$. Note that the choice $N=24$ is usually sufficient, which
is due to the exponential convergence rate ([8, Theorem 4.7]). The next step
is the computation of a singular value decomposition of
$\mathbf{A}_{0,N}=\mathbf{V}\mathbf{\Sigma}\mathbf{W}^{\mathrm{H}}$ with
$\mathbf{V}\in\mathbb{C}^{m\times\ell}$,
$\mathbf{\Sigma}\in\mathbb{C}^{\ell\times\ell}$, and
$\mathbf{W}\in\mathbb{C}^{\ell\times\ell}$. Then, we perform a rank test for
the matrix
$\mathbf{\Sigma}=\mathrm{diag}(\sigma_{1},\sigma_{2},\ldots,\sigma_{\ell})$
for a given tolerance $\epsilon=\text{tol}_{\text{rank}}$ (usually
$\epsilon=10^{-4}$). That is, find $n(\gamma)$ such that
$\sigma_{1}\geq\ldots\geq\sigma_{n(\gamma)}>\epsilon>\sigma_{n(\gamma)+1}\geq\ldots\geq\sigma_{\ell}$.
Define $\mathbf{V}_{0}=(\mathbf{V}_{ij})_{1\leq i\leq m,1\leq j\leq
n(\gamma)}$, $\mathbf{\Sigma}_{0}=(\mathbf{\Sigma}_{ij})_{1\leq i\leq
n(\gamma),1\leq j\leq n(\gamma)}$, and
$\mathbf{W}_{0}=(\mathbf{W}_{ij})_{1\leq i\leq\ell,1\leq j\leq n(\gamma))}$
and compute the $n(\gamma)$ eigenvalues $k_{i}$ and eigenvectors $\vec{s}_{i}$
of the matrix
$\mathbf{B}=\mathbf{V}_{0}^{\mathrm{H}}\mathbf{A}_{1,N}\mathbf{W}_{0}\mathbf{\Sigma}_{0}^{-1}\in\mathbb{C}^{n(\gamma)\times
n(\gamma)}$. The $i$-th non-linear eigenvector $\vec{u}_{i}$ is given by
$\mathbf{V}_{0}\vec{s}_{i}$.

We refer the reader to [8, p. 3849] for more details on the implementation of
this algorithm and the detailed analysis behind it including the proof of
exponential convergence.

###  2.5 Eigenfunction

After we obtain the smallest non-zero eigenvalue $k$ and the corresponding
function $u$ on the boundary from (11), we insert this into (3) to compute the
eigenfunction inside the domain at any point we want. The discretization of
the integrals is done as explained previously. Precisely, we have

$\displaystyle
u(x)=-\mathrm{DL}_{k}^{\Gamma_{1}}u_{1}(x)-\mathrm{DL}_{k}^{\Gamma_{2}}u_{2}(x)\approx\sum_{j=1}^{n_{f}}\sum_{k=1}^{3}\left(\hat{a}_{j,k}u_{1}(\widetilde{v}_{j,k})+\hat{b}_{j,k}u_{2}(\widetilde{v}_{j,k})\right)$

with

$\displaystyle\hat{a}_{j,k}$ $\displaystyle=$
$\displaystyle\int_{\sigma}\partial_{\nu_{1}(\widetilde{m}_{j}(s))}\Phi_{k}(x,\widetilde{m}_{j}(s))\|\partial_{s}\widetilde{m}_{j}(s)\|\widetilde{L}_{k}(s)\,\mathrm{d}s(s)$
$\displaystyle\hat{b}_{j,k}$ $\displaystyle=$
$\displaystyle\int_{\sigma}\partial_{\nu_{2}(\widetilde{m}_{j}(s))}\Phi_{k}(x,\widetilde{m}_{j}(s))\|\partial_{s}\widetilde{m}_{j}(s)\|\widetilde{L}_{k}(s)\,\mathrm{d}s(s)$

for an arbitrary point $x\in D$. In fact, we can find maximal or minimal
values by maximizing or minimizing this function. This is done by the Nelder-
Mead algorithm (in Matlab by the fminsearch function), refer also to [26].

###  2.6 Superconvergence

The convergence of the method is out of the scope of this paper. Standard
convergence results are available for boundary integral equations of the
second kind using boundary element collocation method under suitable
assumptions on the boundary (for example the boundary is at least of class
$C^{2}$) and the boundary condition for the Laplace equation (see [5]).
Quadratic approximations of the boundary and the boundary function yield cubic
convergence (refer to [5] for the Laplace equation and [22] for the Helmholtz
equation). In fact, the convergence results can be improved as shown in [22]
for the three-dimensional case. However, the exact theoretical convergence
rate for the eigenvalue is not known. It is expected that it is at least of
order three, but we see later in the numerical results that it is better than
three (for a sufficiently smooth boundary). Future work in this direction
could be done using the ideas of Steinbach & Unger [38].

Finally, note that a specific choice of $0<\alpha<1/2$ can improve the overall
convergence rate for smooth boundaries (see [22]), but since we are happy with
the cubic convergence rate with the pick $\alpha=(1-\sqrt{3/5})/2$ (a Gauss-
quadrature point within the interval $[0,1]$), we have not investigated this
any further.

##  3 Numerical results

###  3.1 Simply-connected convex domains

First, we check the correctness and the convergence of the underlying method
for the unit circle. It is known that the first non-trivial interior Neumann
eigenvalue is the smallest positive root of $J_{1}^{\prime}$ (the first
derivative of the first kind Bessel function of order one). The root can be
computed to arbitrary precision with Maple with the command

     restart; Digits:=16: fsolve(diff(BesselJ(1,x),x),x=1..2);

It is approximately given by

$1.841\,183\,781\,340\,659\;.$ (12)

For the Beyn algorithm we use the parameters $N=24$, $R=1/2$, $\mu=2$, and
$\ell=10$ for various number of faces $n_{f}$ and number of collocation points
$n_{c}$. With the definition of the absolute error $E_{n_{f}}^{(i)}$ of the
$i$-th eigenvalue approximation, we compute the error of the first non-trivial
eigenvalue $E_{n_{f}}^{(1)}$ of our method compared with (12). Additionally,
we define the estimated error of convergence
$\mathrm{EOC}^{(i)}=\log(E_{n_{f}}^{(i)}/E_{2\cdotp n_{f}}^{(i)})/\log(2)$ of
the $i$-th eigenvalue approximation and compute $\mathrm{EOC}^{(1)}$. As we
can see in Table 1, the first non-trivial Neumann eigenvalue can be made
accurate up to eleven digits.

Table 1:  Absolute error and estimated order of convergence of the first non-
trivial interior Neumann eigenvalue for a unit circle using different number
of faces and collocation points.

  * $n_{f}$ | $n_{c}$ | abs. error $E_{n_{f}}^{(1)}$ | $\mathrm{EOC}^{(1)}$
---|---|---|---
5 | 15 | $5.8503_{-3}$ |
10 | 30 | $4.7818_{-4}$ | 3.6129
20 | 60 | $4.5775_{-5}$ | 3.3849
40 | 120 | $5.0168_{-6}$ | 3.1897
80 | 240 | $5.9173_{-7}$ | 3.0838
160 | 480 | $7.2096_{-8}$ | 3.0369
320 | 960 | $8.9069_{-9}$ | 3.0169
640 | 1920 | $1.1072_{-9}$ | 3.0080
1280 | 3840 | $1.3803_{-10}$ | 3.0039

The estimated order of convergence is at least of order three. Note that we
can go beyond eleven digits accuracy by further increasing $n_{f}$, but it is
not necessary here.

Next, we show the influence of the algebraic multiplicity of the eigenvalue.
The first non-trivial interior Neumann eigenvalue for the unit circle has
algebraic multiplicity two. The same is true for the second non-trivial
interior Neumann eigenvalue. It is obtained by computing the first positive
root of $J^{\prime}_{2}$ approximately given by $3.054\,236\,928\,227\,140$.
The third non-trivial interior Neumann eigenvalue is simple and obtained by
computing the second root of $J_{0}^{\prime}$ given by
$3.831\,705\,970\,207\,512$. Note that the first root of $J_{0}^{\prime}$ is
zero which corresponds to the interior Neumann eigenvalue zero with a
corresponding eigenfunction which is a constant. In Table 2, we show the
absolute error and the estimated order of convergence for the second and third
non-zero interior Neumann eigenvalue for a unit circle where we used different
number of faces and different number of collocation points using the
parameters $N=24$, $R=1/2$, $\mu=3$, and $\ell=10$.

Table 2:  Absolute error and estimated order of convergence of the second and
third non-trivial interior Neumann eigenvalue for a unit circle using
different number of faces and collocation points.

  * $n_{f}$ | $n_{c}$ | abs. error $E_{n_{f}}^{(2)}$ | $\mathrm{EOC}^{(2)}$ | abs. error $E_{n_{f}}^{(3)}$ | $\mathrm{EOC}^{(3)}$
---|---|---|---|---|---
5 | 15 | $1.2817_{-2}$ |  | $1.6335_{-2}$ |
10 | 30 | $1.1543_{-3}$ | 3.4729 | $1.3081_{-3}$ | 3.6424
20 | 60 | $1.2187_{-4}$ | 3.2437 | $1.3133_{-4}$ | 3.3162
40 | 120 | $1.4173_{-5}$ | 3.1041 | $1.5147_{-5}$ | 3.1161
80 | 240 | $1.7199_{-6}$ | 3.0428 | $1.8416_{-6}$ | 3.0401
160 | 480 | $2.1229_{-7}$ | 3.0182 | $2.2788_{-7}$ | 3.0146
320 | 960 | $2.6386_{-8}$ | 3.0082 | $2.8373_{-8}$ | 3.0057
640 | 1920 | $3.2895_{-9}$ | 3.0038 | $3.5406_{-9}$ | 3.0024
1280 | 3840 | $4.1066_{-10}$ | 3.0018 | $4.4225_{-10}$ | 3.0011

Again, we notice the estimated order of convergence of at least order three
and that we can achieve the order and the high accuracy regardless of the
algebraic multiplicity of the eigenvalue.

We define $\square_{\kappa}=[-\kappa,\kappa]\times[-\kappa,\kappa]$ with
$\kappa\in\mathbb{R}_{>0}$. In general, we use a resolution of $100\times 100$
equidistantly distributed points, here within $\square_{1.1}$ and compute for
each point that is located inside the unit circle the value of the
eigenfunction. In Figure 7, we show the eigenfunctions corresponding to the
first three non-trivial interior Neumann eigenvalues as a contour plot with 40
contour lines. We also include the location of the maximum and minimum of the
eigenfunction that corresponds to the first non-trivial interior Neumann
eigenvalue as a red and blue dot, respectively.

(a) First eigenfunction of a unit circle (b) Second eigenfunction of a unit
circle (c) Third eigenfunction of a unit circle

Figure 7:  The first three eigenfunctions corresponding to the first three
non-trivial interior Neumann eigenvalues $1.841\,184$, $3.054\,237$, and
$3.831\,706$ for the unit circle.

We can see that the extreme values for the first non-trivial interior Neumann
eigenfunction of the unit circle are obtained on the boundary as it is
conjectured for simply-connected convex domains. Note that the second
eigenfunction corresponding to the first non-trivial interior Neumann
eigenvalue is a rotated version of the first eigenfunction.

We also show the eigenfunctions including the maximal and minimal value for a
variety of other simply-connected convex domains in Figure 8 such as an
ellipse and two deformed ellipses.

(a) First eigenfunction of an ellipse (b) First eigenfunction of a deformed
ellipse with $\varepsilon=0.1$ (c) First eigenfunction of a deformed ellipse
with $\varepsilon=0.2$

Figure 8:  The eigenfunctions corresponding to the first non-trivial interior
Neumann eigenvalues for the ellipse and two deformed ellipses with
$\varepsilon=0.1$ and $\varepsilon=0.2$ with corresponding non-trivial
interior Neumann eigenvalue $1.544\,422$, $1.849\,064$, and $1.819\,478$.

The boundary of the ellipse is given in parametric form as
$(6\cos(t)/5,\sin(t))^{\top}$ with $t\in[0,2\pi)$. We use the parameters as
before for Beyn’s algorithm except $\mu=1.5$ and consider $\square_{1.3}$. The
first non-trivial interior Neumann eigenvalue is given by $1.544\,422$ which
has algebraic multiplicity one. The parametrization of the deformed ellipse’s
boundary is given by $(0.75\cos(t)+\varepsilon\cos(2t),\sin(t))^{\top}$ with
$t\in[0,2\pi)$, where the parameter $\varepsilon$ is chosen to be $0.1$ and
$0.2$ (see [11, 24] for its first and second use). Using $n_{f}=320$, the
first non-trivial interior Neumann eigenvalue of the deformed ellipses with
$\varepsilon=0.1$ and $\varepsilon=0.2$ are $1.849\,064$ and $1.819\,478$,
respectively. Again, they both have algebraic multiplicity one. It is
generally believed that the hot spots conjecture for general simply-connected
convex domains is true, but a general proof is still open. In all our
numerical results for simply-connected convex domains, we obtain the extrema
on the boundary as one can see in Figure 8.

The same is true when we consider piecewise smooth convex domains such as the
unit square and the equilateral triangle with side length one. The first non-
trivial interior Neumann eigenvalues are known to be $\pi$ and $4\pi/3$ (see
[15]). Their multiplicity is two.

In Table 3 we see that our algorithm works fine with these two piecewise
smooth domains. The estimated order of convergence is better than four and
hence better than for the previously discussed smooth domains. Since we use
quadratic interpolation of the smooth boundary, there is an approximation
error limiting the convergence rate. For the considered piecewise smooth
domains, the boundary is approximated exactly since it is a linear function
thus explaining the better convergence.

Table 3:  Absolute error and estimated order of convergence of the first non-
trivial interior Neumann eigenvalue for a unit square and an equilateral
triangle with side length one using different number of faces and collocation
points.

  * $n_{f}$ | $n_{c}$ | abs. error $E_{n_{f}}^{\square}$ | $\mathrm{EOC}^{\square}$ | $n_{f}$ | $n_{c}$ | abs. error $E_{n_{f}}^{\Delta}$ | $\mathrm{EOC}^{\Delta}$
---|---|---|---|---|---|---|---
4 | 12 | $7.3586_{-3}$ |  | 3 | 9 | $2.7096_{-2}$ |
8 | 24 | $4.4226_{-4}$ | 4.0565 | 6 | 18 | $2.0723_{-3}$ | 3.7088
16 | 48 | $1.9060_{-5}$ | 4.5363 | 12 | 36 | $9.7843_{-5}$ | 4.4046
32 | 96 | $7.6974_{-7}$ | 4.6301 | 24 | 72 | $4.1556_{-6}$ | 4.5573
64 | 192 | $3.0527_{-8}$ | 4.6562 | 48 | 144 | $1.7064_{-7}$ | 4.6060
128 | 384 | $1.2042_{-9}$ | 4.6639 | 96 | 288 | $6.9001_{-8}$ | 4.6282
256 | 768 | $4.7429_{-11}$ | 4.6662 | 192 | 576 | $2.7581_{-10}$ | 4.6449
512 | 1536 | $1.8665_{-12}$ | 4.6674 | 384 | 1152 | $1.0880_{-11}$ | 4.6640

Later, we see that these nice convergence rates depend on the regularity of
the solution at a corner and we obtain worse approximation results. In Figure
9 we show one of the corresponding eigenfunctions for the unit square (refer
also to Figure 1) and the equilateral triangle with side length one including
the location of the maximum and minimum. As we can see, they are located on
the boundary.

(a) First eigenfunction of a unit square (b) First eigenfunction of an
equilateral triangle with side length one

Figure 9:  One of the eigenfunctions corresponding to the first non-trivial
interior Neumann eigenvalues for the unit square and the equilateral triangle
with side length one with corresponding non-trivial interior Neumann
eigenvalue $\pi$ and $4\pi/3$.

###  3.2 Simply-connected non-convex domains

No simply-connected non-convex domain has yet been found (neither
theoretically nor numerically) that fails the hot spots conjecture. Now, we
concentrate on this case. We consider the deformed ellipse from the previous
section with $\varepsilon=0.3$, the peanut-shaped domain and the apple-shaped
domain. The boundaries of the last two domains are given parametrically as
$\sqrt{\cos^{2}(t)+\sin^{2}(t)/4}\left(\cos(t),\sin(t)\right)^{\top}$ and
$\frac{0.5+0.4\cos(t)+0.1\sin(2t)}{1+0.7\cos(t)}\left(\cos(t),\sin(t)\right)^{\top}$
with $t\in[0,2\pi)$, respectively (see [41] for its use). We use the
parameters as before with $\mu=1.5$ for the first two domains and $\mu=3$ for
the apple-shaped domain. Using $n_{f}=320$, the first non-trivial interior
Neumann eigenvalue for the three domains are $1.770\,906$, $1.721\,292$, and
$2.761\,274$, respectively. They are all simple. As we can see in Figure 10,
the maximum and minimum are obtained on the boundary of the domains using
$\square_{1.3}$.

(a) First eigenfunction of a deformed ellipse with $\varepsilon=0.3$ (b) First
eigenfunction of a peanut (c) First eigenfunction of an apple

Figure 10:  The eigenfunctions corresponding to the first non-trivial interior
Neumann eigenvalues for the deformed ellipse with $\varepsilon=0.3$, the
peanut-shaped, and the apple-shaped domain with corresponding non-trivial
interior Neumann eigenvalue $1.770\,906$, $1.721\,292$, and $2.761\,274$.

Interesting domains have been constructed by Kleefeld (refer to [21] for more
details) and extended by Abele and Kleefeld [1] for the purpose of finding new
shape optimizers for certain non-trivial Neumann eigenvalues. The boundaries
of the domains considered in those articles are given by ‘generalized’
equipotentials which are implicit curves. The simplest equipotential is of the
form $\sum_{i=1}^{m}\|x-P_{i}\|^{-1}=c$ where the points
$P_{i}\in\mathbb{R}^{2}$, $i=1,\ldots,m$, the number of points $m$, and the
parameter $c$ are given. All $x\in\mathbb{R}^{2}$ that satisfy the equation
describe the boundary of the domain. We use this idea to construct three non-
symmetric and non-convex simply-connected domains, say $D_{1}$, $D_{2}$, and
$D_{3}$. For the boundary of the domain $D_{1}$, we use the parameter $m=3$
and $c=14/5$. The three points are $(-1,1/2)^{\top}$, $(1,1/3)^{\top}$, and
$(0,4/5)^{\top}$. Using $\mu=1$, we obtain the first non-trivial interior
Neumann eigenvalue $1.051\,055$. For the plot of the eigenfunction, we use
$\square_{1.6}$. The boundary of the second domain is constructed through the
use of $m=4$ and $c=7/2$ and the points are $(-5/4,1/10)^{\top}$,
$(5/4,0)^{\top}$, $(1/10,-1)^{\top}$, and $(0,1)^{\top}$. We obtain the first
non-trivial interior Neumann eigenvalue $1.086\,037$ when using $\mu=1$ and
the eigenfunction using $\square_{1.8}$. The third domain’s boundary is
constructed using $m=5$ and $c=18/5$ and the points are $(-1,1/2)^{\top}$,
$(1,1/2)^{\top}$, $(0,-1)^{\top}$, $(3/2,-1)^{\top}$, and
$(-3/2,-6/5)^{\top}$. Using $\mu=1$, we obtain the first non-trivial interior
Neumann eigenvalue $0.861\,858$. A plot of the corresponding eigenfunction is
shown within $\square_{2.1}$.

All three eigenfunctions including the location of the maximal and minimal
value are shown in Figure 11.

(a) First eigenfunction of $D_{1}$ (b) First eigenfunction of $D_{2}$ (c)
First eigenfunction of $D_{3}$

Figure 11:  The eigenfunctions corresponding to the first non-trivial interior
Neumann eigenvalues for domains $D_{1}$, $D_{2}$, and $D_{3}$ with
corresponding non-trivial interior Neumann eigenvalue $1.051\,055$,
$1.086\,037$, and $0.861\,858$.

As we can see again, the extreme values are obtained on the boundary. The same
is true for the L-shaped domain that is given by
$\llcorner=[0,1]^{2}-[0.5,1]^{2}$. We obtain the first non-trivial interior
Neumann eigenvalue $2.429\,474$ and compare it with the well-known value (see
[14] for the approximation $2\sqrt{1.475\,621\,845}\approx 2.429\,504$). As we
can see in Table 4 the absolute error decreases dramatically since we have
less regularity of the solution at the corner. We are only able to achieve
four digits accuracy with a convergence rate that seems to be $4/3$.

Table 4:  Absolute error and estimated order of convergence of the first non-
trivial interior Neumann eigenvalue for an L-shaped domain using different
number of faces and collocation points.

  * $n_{f}$ | $n_{c}$ | abs. error $E_{n_{f}}^{\llcorner}$ | $\mathrm{EOC}^{\llcorner}$
---|---|---|---
6 | 18 | $7.0712_{-3}$ |
12 | 36 | $3.1014_{-3}$ | 1.1891
24 | 72 | $1.1976_{-3}$ | 1.3728
48 | 144 | $4.7500_{-4}$ | 1.3341
96 | 288 | $1.8860_{-4}$ | 1.3326
192 | 576 | $7.4872_{-5}$ | 1.3329
384 | 1152 | $2.9724_{-5}$ | 1.3327

In Figure 12 we show the corresponding eigenfunction.

(a) First eigenfunction of $\llcorner$ (b) First eigenfunction of
$\llcorner_{2}$ (c) First eigenfunction of $\llcorner_{3}$

Figure 12:  The eigenfunction corresponding to the first non-trivial interior
Neumann eigenvalues for the L-shaped domain and its two variants with
corresponding non-trivial interior Neumann eigenvalue $2.429\,474$,
$2.725\,559$, and $2.207\,276$, respectively.

We also tried two different domains $\llcorner_{2}$ and $\llcorner_{3}$ as
shown in Figure 4 b) and c). Using the same parameters as for the L-shaped
domain $\llcorner$, we obtain the first non-trivial Neumann eigenvalues
$2.725\,559$ and $2.207\,276$, respectively. Interestingly, the approximate
convergence rates are $2.8432$, $2.0744$, $1,6656$, $1.6411$, $1.7140$, and
$2.0115$ for $\llcorner_{2}$ and $1.0716$, $1.1645$, $1.1745$, $1.2202$,
$1.3329$, and $1.6816$ for $\llcorner_{3}$. That is why we concentrate on
smooth boundaries.

In sum, we are not able to construct a simply-connected non-convex domain that
fails the hot spots conjecture. Next, we concentrate on non-simply-connected
domains.

###  3.3 Non-simply-connected domains

Now, we consider an annulus with inner radius $R_{1}=1/2$ and outer radius
$R_{2}=2$. For this domain, we can again compute a reference solution to
arbitrary precision. The non-trivial interior Neumann eigenvalues are obtained
through the roots of

$J_{n}^{\prime}(R_{1}x)Y_{n}^{\prime}(R_{2}x)-J_{n}^{\prime}(R_{2}x)Y_{n}^{\prime}(R_{1}x)=0\,,\quad
n=0,1,2,\ldots$

where $Y_{n}^{\prime}$ denotes the derivative of the second kind Bessel
function of order $n$. (see for example [40, Equation 4.16]). The first two
roots are obtained with the Maple command

     restart; Digits:=16:
     Jp:=unapply(diff(BesselJ(1,x),x),x):
     Yp:=unapply(diff(BesselY(1,x),x),x):
     fsolve(Jp(x/2)*Yp(2*x)-Jp(2*x)*Yp(x/2),x=1);
     Jp:=unapply(diff(BesselJ(2,x),x),x):
     Yp:=unapply(diff(BesselY(2,x),x),x):
     fsolve(Jp(x/2)*Yp(2*x)-Jp(2*x)*Yp(x/2),x=1);

The two smallest roots are approximately given by

$0.822\,252\,688\,623\,884\quad\text{and}\quad 1.504\,647\,782\,189\,479\,,$

respectively. They both have multiplicity two. Again, we show in Table 5 that
our method is able to achieve ten digits accuracy with a cubic convergence
order for the first two non-trivial interior Neumann eigenvalues for an
annulus using the parameter $\mu=6/5$ for various number of faces. Note that
twice the number of faces is needed due to the need of two boundary curves.

Table 5:  Absolute error and estimated order of convergence of the first and
second non-trivial interior Neumann eigenvalue for an annulus with $R_{1}=1/2$
and $R_{2}=2$ using different number of faces and collocation points.

  * $2n_{f}$ | $2n_{c}$ | absolute error $E_{n_{f}}^{(1)}$ | $\mathrm{EOC}^{(1)}$ | absolute error $E_{n_{f}}^{(2)}$ | $\mathrm{EOC}^{(2)}$
---|---|---|---|---|---
10 | 30 | $2.4700_{-3}$ |  | $6.1268_{-3}$ |
20 | 60 | $1.9490_{-4}$ | 3.6637 | $5.4708_{-4}$ | 3.4853
40 | 120 | $1.7741_{-5}$ | 3.4575 | $5.7347_{-5}$ | 3.2540
80 | 240 | $1.8725_{-6}$ | 3.2441 | $6.6443_{-6}$ | 3.1095
160 | 480 | $2.1655_{-7}$ | 3.1122 | $8.0539_{-7}$ | 3.0443
320 | 960 | $2.6159_{-8}$ | 3.0493 | $9.9900_{-8}$ | 3.0111
640 | 1920 | $3.2351_{-9}$ | 3.0154 | $1.2988_{-8}$ | 2.9433

In Figure 13 we show the first three eigenfunctions for $\square_{2.1}$
corresponding to the non-trivial interior Neumann eigenvalues $0.822\,253$,
$1.504\,648$, and $2.096\,773$. To compute the third eigenvalue which has
multiplicity two, we used $\mu=2$. For the first eigenfunction plot, we also
added the extreme values.

(a) First eigenfunction of an annulus (b) Second eigenfunction of an annulus
(c) Third eigenfunction of an annulus

Figure 13:  The first three eigenfunctions corresponding to the first three
non-trivial interior Neumann eigenvalues $0.822\,253$, $1.504\,648$, and
$2.096\,773$ for an annulus.

We can see that the extreme values are again obtained on the boundary of the
annulus. Now, we consider more complex non-simply connected domains. The first
domain $A_{1}$ is given by a unit circle centered at the origin removing an
ellipse centered at $(-0.5,-0.3)^{\top}$ with semi-axis $0.15$ and $0.35$. We
obtain the first non-trivial interior Neumann eigenvalue $1.662\,873$ and the
eigenfunction within $\square_{1.1}$. The second domain $A_{2}$ is given by
$D_{3}$ removing an ellipse centered at $(-0.1,-0.3)^{\top}$ with semi-axis
$0.15$ and $0.35$. We get the first non-trivial interior Neumann eigenvalue
$1.651\,571$ and the eigenfunction within $\square_{1.1}$. The third domain
$A_{3}$ is given by $D_{3}$ removing a 90 degree counter-clockwise rotated
version of $D_{3}$ scaled by $1/2$. We obtain the first non-trivial interior
Neumann eigenvalue $1.171\,590$ and the eigenfunction within $\square_{1.1}$.
We used the parameter $\mu=1.5$ for all three domains. The eigenfunctions for
$A_{1}$, $A_{2}$, and $A_{3}$ including the extreme values are illustrated in
Figure 14.

(a) First eigenfunction of $A_{1}$ (b) First eigenfunction of $A_{2}$ (c)
First eigenfunction of $A_{3}$

Figure 14:  The eigenfunctions corresponding to the first non-trivial interior
Neumann eigenvalues $1.662\,873$, $1.651\,571$, and $1.171\,590$ for the
domains $A_{1}$, $A_{2}$, and $A_{3}$.

As we can see again, the extreme values are obtained on the boundary. We also
tried many other similar geometries with different kinds of holes and obtained
similar results. The hot spots conjecture seems to hold.

Finally, we concentrate on a more complex geometry inspired by the work of
Burdzy [9, Figure 1]. He proved that there exists a bounded planar domain with
one hole that fails the hot spots conjecture. However, the description of the
proposed bounded domain with one hole is very technical and it is difficult to
implement his domain to verify his theoretical result. His domain’s boundary
has many corners which would complicate the explicit construction of the
boundary. Additionally, the domains are very thin. Note that no numerical
result supports his theoretical result. In fact, up-to-date no numerical
results have been shown yet.

We try to close this gap. We construct less complex domains that fail the hot
spots conjecture and show it numerically. Further, we show the exact location
of the extreme values.

The ‘teether’ domain depicted in Figure 15 is constructed as follows:

Figure 15:  A more advance bounded domain with a hole inspired by the work of
Burdzy [9, Figure 1].

Let $E(x,y,a,b,t_{1},t_{2})$ be the ellipse centered at $(x,y)^{\top}$ with
semi-axis $a$ and $b$ constructed for $\phi\in[t_{1},t_{2})$ using the
parametrization $(x+a\cos(\phi),y+b\sin(\phi))^{\top}$. Note that we also
allow $t_{1}>t_{2}$ to guarantee the needed orientation of the curve. The
first half of the outer boundary is given by the pieces
$E(4,4,1,1,-\pi/2,-\pi)$, $E(2,4,1,1,0,\pi)$, $E(0,4,1,1/2,0,-\pi)$,
$E(-2,4,1,1,0,\pi)$, $E(-4,4,1,1,0,-\pi/2)$, and $E(-4,0,3,3,\pi/2,3\pi/2)$.
Rotating this half by $\pi$ yields the second half of the outer boundary. The
first half of the inner boundary is given by the pieces
$E(4,3/2,1,1,\pi/2,\pi)$, $E(2,3/2,1,1,0,-\pi)$, $E(0,3/2,1,1/2,0,\pi)$,
$E(-2,3/2,1,1,0,-\pi)$, $E(-4,3/2,1,1,0,\pi/2)$, and
$E(-4,0,5/2,5/2,\pi/2,3\pi/2)$. Rotating this half by $\pi$ yields the second
half of the inner boundary. Next, the orientation of the inner boundary is
reversed. Finally, all coordinates of the boundary are multiplied with $1/4$.
This yields our ‘teether’ domain $C_{1}$.

We use the parameter $\mu=0.8$, $\ell=20$ and $\square_{1.8}$. We obtain the
first non-trivial interior Neumann eigenvalue $0.370\,708$. The corresponding
eigenfunction including is extreme values is shown in Figure 16. Since the
extreme values lie in very flat plateaus, we additionally show zoomed versions
around the extreme values to better see that they are located inside the
domain.

(a) First eigenfunction of $C_{1}$ (b) Zoom around the maximum (c) Zoom around
the minimum

Figure 16:  The eigenfunction and its zoomed version around the maximum and
minimum corresponding to the first non-trivial interior Neumann eigenvalue
$0.370\,708$ for the teether domain $C_{1}$.

This shows that we are able to show numerically that there exists a bounded
domain with one hole that fails the hot spots conjecture.

Next, we investigate some possible conditions needed to construct an example
that fails the hot spots conjecture. With one bump it was not possible to
obtain the extreme values inside the domain, say $C_{2}$. We use the previous
example, and remove one of the bump and its mirror version in the upper and
lower part of the teether domain. We obtain the first non-trivial interior
Neumann eigenvalue $0.534\,605$ and the corresponding eigenfunction within
$\square_{1.6}$.

(a) First eigenfunction of $C_{2}$ (b) Zoom around the maximum

Figure 17:  The eigenfunction and its zoomed version around the maximum
corresponding to the first non-trivial interior Neumann eigenvalue
$0.534\,605$ for the domain $C_{2}$.

As we can see again in Figure 17, the values inside the bump area are very
close to each other. The extreme values are attained on the boundary.

Right now, the proposed domain $C_{1}$ has two lines of symmetry. Now, we
break the symmetry and show that we are still able to obtain a counter-example
of the hot spots conjecture. We add a small amount of $0.0250$ to the semi
axis of the ellipse that describes the upper left bump and obtain the domain
$C_{3}$. The first non-trivial interior Neumann eigenvalue is $0.367\,496$.

(a) First eigenfunction of $C_{3}$ (b) First eigenfunction of
$\widetilde{C}_{3}$

Figure 18:  The eigenfunction corresponding to the first non-trivial interior
Neumann eigenvalue $0.367\,496$ and $0.370\,054$ for the domains $C_{3}$ and
$\widetilde{C}_{3}$.

As we can see in Figure 18, the location of the hot spots of the eigenfunction
within $\square_{1.8}$ are slightly changed, too, but they remain inside of
$C_{3}$. Increasing the value from $0.0250$ to $0.1255$ will shift the maximal
value from inside the domain $\widetilde{C}_{3}$ to the boundary while the
minimum stays inside the domain. We obtain the eigenvalue $0.370\,054$.

Now, we show what happens if we make the gap between $E(0,4,1,1/2,0,-\pi)$ and
$E(0,3/2,1,1/2,0,\pi)$ and its mirror counterpart smaller. We use
$E(0,4,1,1,0,-\pi)$ and $E(0,3/2,1,1,0,\pi)$ instead to obtain $C_{4}$. We
obtain the first non-trivial interior Neumann eigenvalue $0.384\,715$ and its
corresponding eigenfunction within $\square_{1.8}$ shown in Figure 19.

(a) First eigenfunction of $C_{4}$ Figure 19:  The eigenfunction corresponding
to the first non-trivial interior Neumann eigenvalue $0.384\,715$ for the
domains $C_{4}$.

Next, we show what happens if we make the bumps of the domain $C_{4}$ smaller.
We construct the domain as before except that we introduce a new parameter
$\delta>0$. The first half of the outer boundary is given by the pieces
$E(4,3+\delta,1,\delta,-\pi/2,-\pi)$, $E(2,3+\delta,1,\delta,0,\pi)$,
$E(0,3+\delta,1,\delta,0,-\pi)$, $E(-2,3+\delta,1,\delta,0,\pi)$,
$E(-4,3+\delta,1,\delta,0,-\pi/2)$, and $E(-4,0,3,3,\pi/2,3\pi/2)$. Rotating
this half by $\pi$ yields the second half of the outer boundary. The first
half of the inner boundary is given by the pieces
$E(4,5/2-\delta,1,\delta,\pi/2,\pi)$, $E(2,5/2-\delta,1,\delta,0,-\pi)$,
$E(0,5/2-\delta,1,\delta,0,\pi)$, $E(-2,5/2-\delta,1,\delta,0,-\pi)$,
$E(-4,5/2-\delta,1,\delta,0,\pi/2)$, and $E(-4,0,5/2,5/2,\pi/2,3\pi/2)$.
Rotating this half by $\pi$ yields the second half of the inner boundary.
Next, the orientation of the inner boundary is reversed. Finally, all
coordinates of the boundary are multiplied with $1/4$. Note that $\delta=1$
yields the domain $C_{4}$.

We obtain the following results within $\square_{1.8}$ shown in Figure 20 for
$\delta=1/4$, $\delta=1/10$, and $\delta=1/20$.

(a) First eigenfunction of $\tilde{C}_{4}$ (b) First eigenfunction of
$\hat{C}_{4}$ (c) First eigenfunction of $\bar{C}_{4}$

Figure 20:  The eigenfunction corresponding to the first non-trivial interior
Neumann eigenvalue $0.578\,402$, $0.668\,373$, and $0.708\,633$ for the
domains $\tilde{C}_{4}$, $\hat{C}_{4}$ and $\bar{C}_{4}$, respectively.

Surprisingly, the extreme values stay inside of the domain $\tilde{C}_{4}$,
$\hat{C}_{4}$, and $\bar{C}_{4}$. Also note that the first non-trivial
interior Neumann eigenvalue changes drastically.

Interestingly, we can also remove the bumps in the lower part of the teether
domain $C_{1}$ and still obtain the extreme values inside the new domain
$C_{5}$. The result is shown in Figure 21 within $\square_{1.8}$. The
corresponding non-trivial interior Neumann eigenvalue is $0.563\,329$.

(a) First eigenfunction of $C_{5}$ Figure 21:  The eigenfunction corresponding
to the first non-trivial interior Neumann eigenvalue $0.563\,329$ for the
domain $C_{5}$.

The idea to also remove the bumps on the upper part yields the following
results for the new ‘stadium’ domain $S$ as shown in Figure 22 within
$\square_{1.8}$.

(a) First eigenfunction of $S$ (b) Second eigenfunction of $S$

Figure 22:  The eigenfunctions corresponding to the first and second non-
trivial interior Neumann eigenvalue $0.755\,416$ and $0.756\,082$ for the
domain $S$.

Unfortunately, the extreme values are now on the boundary of $S$. The first
non-trivial interior Neumann eigenvalues is $0.755\,416$. The second non-
trivial eigenvalue is very close (but distinct) $0.756\,082$. We can see that
the limiting process of making the bumps of $C_{4}$ smaller (see also
$\tilde{C}_{4}$ and $\hat{C}_{4}$) yields the eigenfunction corresponding to
the second non-trivial interior Neumann eigenfunction for the domain $S$. This
is very unexpected.

Further, it seems that the heat flow out of the narrow part of the ‘pipes’ has
to be almost without inclination. We use the parametrization
$r_{i}(t)\cdotp(\sin(t),\cos(t))^{\top}$ with $r_{i}(t)=a_{i}(1+1/2\sin(\omega
t\cdotp\mathbb{I}_{[(\omega/2-1)\pi/\omega,(\omega/2+1)\pi/\omega]\cup[(3\omega/2-1)\pi/\omega,(3\omega/2+1)\pi/\omega]}))$,
$t\in[0,2\pi)$ with $a_{1}=1$ for the outer boundary and $a_{2}=0.9$ for the
inner boundary. Here $\mathbb{I}$ denotes the indicator function. We use
$\omega=4$ and $\omega=8$ to construct the domains $F_{1}$ and $F_{2}$. We
obtain the first non-trivial interior Neumann eigenvalue $1.592\,787$ and
$1.717\,098$ and the eigenfunction within $\square_{1.6}$. As we can see in
Figure 23 the maximum and minimum value are obtained on the boundary. In fact,
we have two maxima and two minima.

(a) First eigenfunction of $F_{1}$ (b) First eigenfunction of $F_{2}$

Figure 23:  The eigenfunction corresponding to the first non-trivial interior
Neumann eigenvalue $1.592\,787$ and $1.717\,098$ for the domains $F_{1}$ and
$F_{2}$, respectively.

Next, we show in Table 6 for the different examples that show a failure of the
hot spots conjecture the following information: The domain under
consideration, the location of the global maximum and minimum inside the
domain, and the ratio $\aleph_{\textit{max}}$ and $\aleph_{\textit{min}}$
which are defined as the ratio of the maximum inside the domain divided by the
the maximum on the boundary and likewise for the minimum. All ratios will be
larger than one, but we are interested in how close they are to one. Recall
that the hot spots conjecture fails if we have $\aleph_{\textit{max}}>1$
and/or $\aleph_{\textit{min}}>1$. The maximum and minimum are calculated
through the Matlab minimization routine fminsearch by passing the function
given in (2) to a tolerance of $10^{-10}$ for the step size and the difference
of a function evaluation. For the first two domains, we use as starting value
$(0,0.6)^{\top}$ and $(0,-0.6)^{\top}$, respectively. For the next three
domains, we use $(0,0.69)^{\top}$ and $(0,-0.69)^{\top}$ and for the last
domain, we use $(0,0.6875)^{\top}$ and $(0,-0.6875)^{\top}$.

Table 6:  Location of the maximum and minimum inside the domain along with the
ratios $\aleph_{\textit{max}}>1$ and/or $\aleph_{\textit{min}}>1$ for various
domains $D$ that fail the hot spots conjecture.

  * $D$ | location max | location min | $\aleph_{\textit{max}}$ | $\aleph_{\textit{min}}$
---|---|---|---|---
$C_{1}$ | $(-3.805_{-8},6.877_{-1})^{\top}$ | $(3.299_{-8},-6.877_{-1})^{\top}$ | $1+1.221_{-4}$ | $1+1.221_{-4}$
$C_{4}$ | $(8.126_{-8},6.875_{-1})^{\top}$ | $(-3.893_{-8},-6.875_{-1})^{\top}$ | $1+1.196_{-5}$ | $1+1.196_{-5}$
$\tilde{C}_{4}$ | $(-2.998_{-8},6.875_{-1})^{\top}$ | $(-2.028_{-8},-6.875_{-1})^{\top}$ | $1+7.130_{-6}$ | $1+7.118_{-6}$
$\hat{C}_{4}$ | $(3.530_{-8},6.875_{-1})^{\top}$ | $(2.824_{-8},-6.875_{-1})^{\top}$ | $1+3.808_{-6}$ | $1+3.802_{-6}$
$\bar{C}_{4}$ | $(-4.580_{-8},6.875_{-1})^{\top}$ | $(-4.787_{-8},-6.875_{-1})^{\top}$ | $1+2.137_{-6}$ | $1+2.138_{-6}$
$C_{5}$ | $(-1.968_{-7},6.877_{-1})^{\top}$ | ————— | $1+1.438_{-3}$ | —————

As we can see the maximum and minimum are approximately located on the
$y$-axis each centered between the bumps for the first four domains.
Interestingly, we obtain for the ratios $1+\epsilon$ with very small
$\epsilon>0$. The parameter $\epsilon$ decreases as the shape gets closer to
the ’stadium’ domain. Hence, we might conjecture that there exists a domain
with one hole where we have $\aleph_{\textit{max}}=1+\epsilon$ and
$\aleph_{\textit{min}}=1+\epsilon$ with $\epsilon\geq\epsilon_{0}>0$ given
$\epsilon_{0}$ arbitrarily small. For the last domain, we obtain the largest
value for $\epsilon$. Precisely, we obtain $\epsilon_{0}>1+10^{-3}$. However,
the location of the minimum might be on the boundary (or on the complete part
of the $y$-axis).

Finally, we also show a different easy to construct domain with one hole that
fails the hot spots conjecture. The outer boundary is given by straight lines
connecting the $21$ points $(0,1)^{\top}$, $(3,1)^{\top}$, $(3,0)^{\top}$,
$(4,0)^{\top}$, $(4,1)^{\top}$, $(6,1)^{\top}$, $(6,0)^{\top}$,
$(7,0)^{\top}$, $(7,1)^{\top}$, $(10,1)^{\top}$, $(10,6)^{\top}$,
$(7,6)^{\top}$, $(7,7)^{\top}$, $(6,7)^{\top}$, $(6,6)^{\top}$,
$(4,6)^{\top}$, $(4,7)^{\top}$, $(3,7)^{\top}$, $(3,6)^{\top}$,
$(0,6)^{\top}$, and $(0,1)^{\top}$. The inner boundary is given by the
straight lines connecting the $21$ points $(1,2)^{\top}$, $(3,2)^{\top}$,
$(3,3)^{\top}$, $(4,3)^{\top}$, $(4,2)^{\top}$, $(6,2)^{\top}$,
$(6,3)^{\top}$, $(7,3)^{\top}$, $(7,2)^{\top}$, $(9,2)^{\top}$,
$(9,5)^{\top}$, $(7,5)^{\top}$, $(7,4)^{\top}$, $(6,4)^{\top}$,
$(6,5)^{\top}$, $(4,5)^{\top}$, $(4,4)^{\top}$, $(3,4)^{\top}$,
$(3,5)^{\top}$, $(1,5)^{\top}$, and $(1,2)^{\top}$ and then stored in reversed
order. The resulting coordinates are first shifted by $-(5,7/2)^{\top}$ and
then scaled by $1/2$ to center the ‘brick’ domain $B$ with respect to the
origin. Using $600$ collocation points with the parameters $N=24$, $R=1/2$,
$\mu=7/10$, and $\ell=20$ yields the eigenfunction shown in Figure 24 within
$\square_{3}$ using a resolution of $100\times 100$.

(a) First eigenfunction of $B$ Figure 24:  The eigenfunction corresponding to
the first non-trivial interior Neumann eigenvalue $0.411\,448$ for the brick
domains $B$.

As we can see, the maximal and minimal value are attained within the domain.
Hence, we constructed another domain with one hole that fails the hot spots
conjecture.

###  3.4 Domains with more than one hole

Finally, we show without further discussion that we are also able to construct
examples with more than one hole where the hot spots conjecture fails to hold.
Using the teether domain $C_{1}$ and removing a circle centered at
$(1/2,11/16)^{\top}$ with radius $\mathfrak{R}$, yields a domain with two
holes, say $C_{1,\mathfrak{R}}$. Using $\mathfrak{R}=0.1$,
$\mathfrak{R}=0.15$, and $\mathfrak{R}=0.2$, gives the results shown in Figure
25 where we used the same set of parameters as for the results for $C_{1}$.

(a) First eigenfunction of $C_{1,0.1}$ (b) First eigenfunction of $C_{1,0.15}$
(c) First eigenfunction of $C_{1,0.2}$

Figure 25:  The eigenfunction corresponding to the first non-trivial interior
Neumann eigenvalue $0.372\,580$, $0.374\,917$, and $0.378\,131$ for the
domains $C_{1,0.1}$, $C_{1,0.15}$ and $C_{1,0.2}$, respectively.

As we can see, the maximum and minimum still remain inside the domains
$C_{1,0.1}$ and $C_{1,0.15}$. The location of the maximum is slightly shifted
to the right whereas the minimum is slightly shifted to the left for the first
two cases. If the radius of the removed circle is large, then the maximum goes
to the boundary whereas the minimum stays inside the domain $C_{1,0.2}$ as
shown in the last contour plot of Figure 25. Next, we construct domains with
three holes. Therefore, we use the previous domain $C_{1,\mathfrak{R}}$ and
mirror the circular hole at the $y$-axis. This yields the domain
$\widetilde{C}_{1,\mathfrak{R}}$. Using the same parameters as before yields
the results shown in Figure 26.

(a) First eigenfunction of $\widetilde{C}_{1,0.1}$ (b) First eigenfunction of
$\widetilde{C}_{1,0.15}$ (c) First eigenfunction of $\widetilde{C}_{1,0.2}$

Figure 26:  The eigenfunction corresponding to the first non-trivial interior
Neumann eigenvalue $0.374\,552$, $0.379\,670$, and $0.387\,461$ for the
domains $\widetilde{C}_{1,0.1}$, $\widetilde{C}_{1,0.15}$ and
$\widetilde{C}_{1,0.2}$, respectively.

As we can see now, the maximum and minimum remain inside all the three domains
$\widetilde{C}_{1,0.1}$, $\widetilde{C}_{1,0.15}$, and
$\widetilde{C}_{1,0.2}$. Next, we consider domains with four holes. Therefore,
we use the domains $C_{1,0.1}$, $C_{1,0.15}$ and $C_{1,0.2}$ and mirror the
upper right circular hole at the $x$-axis. This yields the domains
$C_{1,0.1}^{\mathrm{mir}}$, $C_{1,0.15}^{\mathrm{mir}}$ and
$C_{1,0.2}^{\mathrm{mir}}$, respectively. The results are presented in Figure
27.

(a) First eigenfunction of $C_{1,0.1}^{\mathrm{mir}}$ (b) First eigenfunction
of $C_{1,0.15}^{\mathrm{mir}}$ (c) First eigenfunction of
$C_{1,0.2}^{\mathrm{mir}}$

Figure 27:  The eigenfunction corresponding to the first non-trivial interior
Neumann eigenvalue $0.376\,406$, $0.383\,781$, and $0.394\,522$ for the
domains $C_{1,0.1}^{\mathrm{mir}}$, $C_{1,0.15}^{\mathrm{mir}}$ and
$C_{1,0.2}^{\mathrm{mir}}$, respectively.

As expected, we obtain the minimal and maximal value inside the domain for the
two domains $\widetilde{C}_{1,0.1}^{\mathrm{mir}}$ and
$\widetilde{C}_{1,0.15}^{\mathrm{mir}}$ whereas the maximum is on the boundary
for the domain $\widetilde{C}_{1,0.2}^{\mathrm{mir}}$. The domains with five
holes $\widetilde{C}_{1,0.1}^{\mathrm{mir}}$,
$\widetilde{C}_{1,0.15}^{\mathrm{mir}}$ and
$\widetilde{C}_{1,0.2}^{\mathrm{mir}}$ are constructed by mirroring the
domains $\widetilde{C}_{1,0.1}$, $\widetilde{C}_{1,0.15}$ and
$\widetilde{C}_{1,0.2}$ at the $x$-axis. The results are shown in Figure 28.

(a) First eigenfunction of $\widetilde{C}_{1,0.1}^{\mathrm{mir}}$ (b) First
eigenfunction of $\widetilde{C}_{1,0.15}^{\mathrm{mir}}$ (c) First
eigenfunction of $\widetilde{C}_{1,0.2}^{\mathrm{mir}}$

Figure 28:  The eigenfunction corresponding to the first non-trivial interior
Neumann eigenvalue $0.378\,360$, $0.388\,438$, and $0.403\,504$ for the
domains $\widetilde{C}_{1,0.1}^{\mathrm{mir}}$,
$\widetilde{C}_{1,0.15}^{\mathrm{mir}}$ and
$\widetilde{C}_{1,0.2}^{\mathrm{mir}}$, respectively.

As we can see, the maximal and minimal values are attained inside all the
considered domains with five holes.

The extension for the construction of domains having more than five holes
which do not satisfy the hot spots conjecture is now straightforward.

##  4 Summary and outlook

In this paper, a detailed description is given on how to compute the first
non-trivial eigenvalue and its corresponding eigenfunction for the Laplace
equation with Neumann boundary condition for a given domain with one hole. The
problem is reformulated as a non-linear eigenvalue problem involving boundary
integral equations thus reducing a two-dimensional problem to a one-
dimensional problem. Due to superconvergence we are able to achieve highly
accurate approximations both for the eigenvalue and the eigenfunction. With
this method at hand, we can compute the eigenvalue and eigenfunction for
several different constructed domains. This gives the possibility to find
domains with one hole failing the hot spots conjecture and investigate the
influence of varying the domain. Some interesting observation can be made such
as that the ratio between the maximal/minimal value inside the domain and the
maximal/minimal value on the boundary can be $1+10^{-3}$. The Matlab codes
including the produced data are available at github

https://github.com/kleefeld80/hotspots

and researchers can run it on their own constructed domains and reproduce the
numerical results within this article. This might give new ideas whether one
can find assumptions in order to prove or disprove the hot spots conjecture.
The extension for domains with more than one hole is straightforward. For the
sake of completeness they are given at the end of the numerical results
section for domains with up to five holes, but without detailed discussion.

It would be interesting to check whether it is possible to construct three-
dimensional domains with one hole that fail the hot spots conjecture, too. The
software for a domain without a hole would already be available and only needs
to be extended (see [22, 20]). The consideration of other partial differential
equations in two or three dimensions whose fundamental solution is known
together with Neumann boundary condition could be numerically investigated as
well.

The author thanks Prof. Stefan Steinerberger from the University of
Washington, Seattle (USA) for the fruitful discussions during the preparation
of the manuscript.

## References

## References

  * [1] D. Abele and A. Kleefeld.  New numerical results for the optimization of Neumann eigenvalues.  In C. Constanda, editor, Computational and Analytic Methods in Science and Engineering, pages 1–20. Birkhäuser, 2020.
  * [2] E. O. Asante-Asamani, A. Kleefeld, and B. A. Wade.  A second-order exponential time differencing scheme for non-linear reaction-diffusion systems with dimensional splitting.  J. Comput. Phys., 415:109490, 2020.
  * [3] R. Atar.  Invariant wedges for a two-point reflecting Brownian motion and the “hot spots” problem.  Electronic Journal of Probability, 6(18):1–19, 2001.
  * [4] R. Atar and K. Burdzy.  On Neumann eigenfunctions in lip domains.  Journal of the American Mathematical Society, 17(2):243–265, 2004\.
  * [5] K. E. Atkinson.  The Numerical Solution of Integral Equations of the Second Kind.  Cambridge University Press, 1997.
  * [6] R. Bañuelos and K. Burdzy.  On the “hot spots” conjecture of J. Rauch.  Journal of Functional Analysis, 164:1–33, 1999.
  * [7] R. F. Bass and K. Burdzy.  Fiber Brownian motion and the “hot spots” problem.  Duke Mathematical Journal, 105(1):25–58, 2000.
  * [8] W.-J. Beyn.  An integral method for solving nonlinear eigenvalue problems.  Linear Algebra and its Applications, 436:3839–3863, 2012.
  * [9] K. Burdzy.  The hot spots problem in planar domains with one hole.  Duke Mathematical Journal, 129(3):481–502, 2005.
  * [10] K. Burdzy and W. Werner.  A counterexample to the “hot spots” conjecture.  Annals of Mathematics, 149(1):309–317, 1999.
  * [11] F. Cakoni and R. Kress.  A boundary integral equation method for the transmission eigenvalue problem.  Applicable Analysis, 96(1):23–38, 2017.
  * [12] D. Colton and R. Kress.  Inverse acoustic and electromagnetic scattering theory.  Springer, 3rd edition, 2013.
  * [13] P. Freitas.  Closed nodal lines and interior hot spots of the second eigenfunction of the Laplacian on surfaces.  Indiana University Mathematics Journal, 51(2):305–316, 2002.
  * [14] A. Gilette, C. Gross, and K. Plackowski.  Numerical studies of serendipity amd tensor product elements for eigenvalue problems.  Involve: A Journal of Mathematics, 11(4):661–678, 2018.
  * [15] D. S. Grebenkov and B.-T. Nguyen.  Geometrical structure of Laplacian eigenfunctions.  SIAM Review, 55(4):601–667, 2013.
  * [16] R. Hempel, L. A. Seco, and B. Simon.  The essential spectrum of Neumann Laplacians on some bounded singular domains.  Journal of Functional Analysis, 102(2):448–483, 1991.
  * [17] D. Jerison and N. Nadirashvili.  The “hot spots” conjecture for domains with two axes of symmetry.  Journal of the American Mathematical Society, 13(4):741–772, 2000\.
  * [18] C. Judge and S. Mondal.  Euclidean triangles have no hot spots.  Annals of Mathematics, 191(1):167–211, 2020.
  * [19] B. Kawohl.  Rearrangements and Convexity of Level Sets in PDE.  Lecture Notes in Mathematics. Springer, 1985.
  * [20] A. Kleefeld.  Numerical methods for acoustic and electromagnetic scattering: Transmission boundary-value problems, interior transmission eigenvalues, and the factorization method.  Habilitation thesis, Brandenburg University of Technology Cottbus - Senftenberg, Cottbus, 2015.
  * [21] A. Kleefeld.  Shape optimization for interior Neumann and transmission eigenvalues.  In C. Constanda and P. Harris, editors, Integral Methods in Science and Engineering, pages 185–196. Springer, 2019.
  * [22] A. Kleefeld and T.-C. Lin.  Boundary element collocation method for solving the exterior Neumann problem for Helmholtz’s equation in three dimensions.  Electronic Transactions on Numerical Analysis, 39:113–143, 2012\.
  * [23] A. Kleefeld and T.-C. Lin.  A global Galerkin method for solving the exterior Neumann problem for the Helmholtz equation using Panich’s integral equation approach.  SIAM Journal on Scientific Compututing, 35(3):A1709–A1735, 2013\.
  * [24] A. Kleefeld and L. Pieronek.  The method of fundamental solutions for computing acoustic interior transmission eigenvalues.  Inverse Problems, 34(3):035007, 2018.
  * [25] D. Krejčiřík and M. Tušek.  Location of hot spots in thin curved strips.  Journal of Differential Equations, 266(6):2953–2977, 2019.
  * [26] J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright.  Convergence properties of the Nelder–Mead simplex method in low dimensions.  SIAM Journal on Optimization, 9(1):112–147, 1998.
  * [27] R. R. Lederman and S. Steinerberger.  Extreme values of the Fiedler vector on trees.  arXiv 1912.08327, 2019.
  * [28] W. McLean.  Strongly Elliptic Systems and Boundary Integral Operators.  Cambridge University Press, 2000.
  * [29] Y. Miyamoto.  The “hot spots” conjecture for a certain class of planar convex domains.  Journal of Mathematical Physics, 50(10):103530, 2009.
  * [30] Y. Miyamoto.  A planar convex domain with many isolated “hot spots” on the boundary.  Japan Journal of Industrial and Applied Mathematics, 30:145–164, 2013.
  * [31] M. N. Pascu.  Scaling coupling of reflecting Brownian motions and the hot spots problem.  Transactions of the American Mathematical Society, 354(11):4681–4702, 2002.
  * [32] J. Rauch.  Lecture #1. Five problems: An introduction to the qualitative theory of partial differential equations.  In J. Goldstein, editor, Partial differential equations and related topics, volume 446 of Lecture Notes in Mathematics, pages 355–369. Springer, 1974.
  * [33] M. Reed and B. Simon.  Methods of Modern Mathematical Physics: Vol.: 4. : Analysis of Operators.  Academic Press, 1978.
  * [34] S. Sauter and C. Schwab.  Boundary Element Methods, volume 39 of Computational Mathematics.  Springer, 2011.
  * [35] A. F. Seybert, B. Soenarko, F. J. Rizzo, and D. J. Shippy.  An advance computational method for radiation and scattering of acoustic waves in three dimensions.  Journal of the Acoustical Society of America, 77(2):362–368, 1985\.
  * [36] L. F. Shampine.  Vectorized adaptive quadrature in MATLAB.  Journal of Computational and Applied Mathematics, 211(2):131–140, 2008.
  * [37] B. Siudeja.  Hot spots conjecture for a class of acute triangles.  Mathematische Zeitschrift, 208:783–806, 2015.
  * [38] O. Steinbach and G. Unger.  Convergence analysis of a Galerkin boundary element method for the Dirichlet Laplacian eigenvalue problem.  SIAM Journal on Numerical Analysis, 50(2):710–728, 2012.
  * [39] S. Steinerberger.  Hot spots in convex domains are in the tips (up to an inradius).  Communications in Partial Differential Equations, 45(6):641–654, 2020.
  * [40] C. C. Tsai, D. L. Young, C. W. Chen, and C. M. Fan.  The method of fundamental solutions for eigenproblems in domains with and without interior holes.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 462(2069):1443–1466, 2006.
  * [41] J. Yang, B. Zhang, and H. Zhang.  The factorization method for reconstructing a penetrable obstacle with unknown buried objects.  SIAM Journal on Applied Mathematics, 73(2):617–635, 2013.