# A counterexample to Payne’s nodal line conjecture with few holes

Joel Dahne  Javier Gómez-Serrano  Kimberly Hou

###### Abstract

Payne conjectured in 1967 that the nodal line of the second Dirichlet
eigenfunction must touch the boundary of the domain. In their 1997
breakthrough paper, Hoffmann-Ostenhof, Hoffmann-Ostenhof and Nadirashvili
proved this to be false by constructing a counterexample in the plane with
many holes and raised the question of the minimum number of holes a
counterexample can have. In this paper we prove it is at most 6.

##  1 Introduction

Let $\Omega$ be a bounded planar domain, and let $\lambda_{i}$ be the
eigenvalues of the Dirichlet Laplacian, namely the numbers
$0<\lambda_{1}<\lambda_{2}\leq\lambda_{3}\leq\ldots$ that satisfy

$\displaystyle-\Delta u_{k}$ $\displaystyle=\lambda_{k}u_{k}\text{ in }\Omega$
$\displaystyle u_{k}$ $\displaystyle=0\text{ on }\partial\Omega.$ (1)

The so-called “nodal line conjecture” from 1967 by Payne [44, Conjecture 5,
p.467], [45] says that the nodal line (the zero level set) of $u_{2}$ on a
bounded domain in $\mathbb{R}^{2}$ must touch the boundary. This statement was
later extended by Yau [56] to the higher dimensional case. Contrary to its
apparent simplicity, the conjecture is still open in the general case, though
a few results (both positive and negative) and extensions have been done. If
the domain is convex, Melas [41] (in the case of $C^{\infty}$ boundary) and
Alessandrini [1] (in the general case) proved it in the positive. Jerison [30]
proved the conjecture for long thin convex sets, and Jerison [31],
Grieser–Jerison [27] and Beck–Canzani–Marzuola [6] gave more information on
the location of the nodal line. Under various additional symmetry and/or
convexity assumptions, the conjecture has been proved by Payne [45], Lin [37],
Pütter [48], Damascelli [14] and Yang–Guo [55]. See also [26].

In 1997, Hoffmann-Ostenhof–Hoffmann-Ostenhof–Nadirashvili [29] constructed a
counterexample of a planar, bounded, non-simply connected domain for which the
nodal line is closed and does not touch the boundary. Their construction was
extended by Fournais [16] and later by Kennedy [34] for the higher dimensional
case. Freitas–Krejčiřík [18] constructed an unbounded counterexample. Kennedy
[35] proved similar results (both positive and negative) for a toy Neumann
analogue under the assumption of central symmetry of the domain. In both [29,
16] the main domain construction is based on “carving” a sufficiently large
number $N$ of holes on a symmetric domain so that the nodal set gets
disconnected from the boundary, but no quantitative estimates are provided on
$N$ (Kennedy’s construction [34] of a simply connected domain is valid in
dimensions 3 or higher). Specifically, in [29] the number of holes is claimed
to be delicate to bound and in [16] of the order of $10^{9}$. Hoffmann-
Ostenhof, Hoffmann-Ostenhof and Nadirashvili raised the following question
circa 25 years ago:

######  Question 1.1 ([29], Remark 3).

Clearly, an interesting question is whether there exists a simply connected
domain for which the second eigenfuntion has a closed nodal line. We do not
believe this. So a more general question is: What is the smallest $N_{0}$ such
that there exists a domain with $N_{0}$ boundary components whose second
eigenfunction has a nodal line that does not hit the boundary?

In this paper we provide a partial answer to Question 1.1, giving an upper
bound of 7 by constructing an example of a domain with a closed nodal line
(see Figure 1(a) for an illustration of the domain and Figure 1(b) for the
nodal line). Our main result is the following theorem:

######  Theorem 1.2.

There exists a planar domain with 6 holes for which the nodal line of $u_{2}$
is closed.

######  Remark 1.3.

We make no claim that $N_{0}=7$ is optimal. Possibly with a more thorough
search, more computational time and better estimates one can get to $N_{0}=5$
or lower.

We summarize the main steps of the proof. The underlying idea is that if we
can find a suitable candidate domain and we can find an approximate $u_{2}$
(in the sense that it satisfies equation (1) up to a small error) which has a
closed nodal line, then there has to exist an exact solution nearby also
having a closed nodal line, and no other spurious nodal domains can appear
since $u_{2}$ can have at most two nodal domains. The difficult part is to
derive effective (as opposed to asymptotic, either with respect to the number
of holes or without explicit constants) stability estimates that are good
enough, as well as to find good candidates since the problem seems to be quite
unstable in that regard. Moreover, since the candidate domains are far from
any explicit domains for which the spectrum and the nodal lines are known,
perturbative methods fail. Furthermore, the errors need to be very small due
to the nodal line being very close to the boundary. On top of that, the
eigenvalue $\lambda_{2}$ is very close to $\lambda_{3}$ and $\lambda_{4}$
which makes difficult to distinguish the eigenfunction associated with the
former from the latter two. To overcome these difficulties parts of the proof
will be computer-assisted. Nevertheless, one has to proceed very carefully and
derive extremely tight bounds to make all the estimates work.

###  1.1 The Method of Particular solutions and Validation of Eigenvalues

Perturbative analysis of the eigenpairs is a classical problem, and there is
an extensive literature [33, 36, 7, 17, 42, 51, 4] showing the existence of an
eigenvalue close to an approximate one. The theme of these results is that if
one can find $(\lambda_{\text{app}},u_{\text{app}})$ satisfying the equation
up to an error bounded by $\delta$, then one can show that there is a true
eigenpair $(\lambda,u)$ at a distance $C\delta^{k}$. The main drawback is that
these methods can not tell the position of the eigenvalue within the spectrum.

To solve this issue, Plum [46] proposed a homotopy method linking the
eigenvalues of the domain of the problem with another, known domain (the base
problem). See also the intermediate method [54, 20, 5] for another example of
connecting the problem to a known domain. Our domain is very far from any
known domain (due to the holes), so we have opted for a more direct approach.
Using domain monotonicity with inclusion also did not yield good enough
(lower) bounds for our purpose. Using the Finite Element framework, Liu and
Liu–Oishi [40], [38] have managed to give explicit, rigorous computable lower
bounds of the spectrum in terms of solutions of a (big) finite linear system
(see also [11] for similar bounds, [57] for the case of the Steklov problem
and [39] for a more general setting). Nonetheless, the aforementioned bounds
are not good enough to obtain Theorem 1.2 as the eigenvalues are tightly
clustered and that would require a mesh so refined we would not be able to
handle it on a computer within a reasonable time.

In this paper we will combine the two families as in [23]. The first pass will
separate the first 4 eigenvalues from the rest of the spectrum (using the
method of [38]). The enclosures and the scale are coarse at this point. The
second pass will find 4 approximate eigenpairs below the threshold and use the
finer stability methods to separate $\lambda_{2}$ from the others. There is a
big technical difficulty since $\lambda_{3}$ and $\lambda_{4}$ presumably
correspond to a double eigenvalue, and we have to handle the circumstance that
the span is two-dimensional in that case.

In order to find accurate approximations of the eigenvalues and eigenfunctions
we will use the Method of Particular Solutions (MPS). This method was
introduced by Fox, Henrici and Moler [17] and has been later adapted by many
authors (see [2, 49, 21, 15, 8] as a sample, and the thorough review [9]). The
main idea is to consider a set of functions that solve the eigenvalue problem
without boundary conditions as a basis, and writing the solution of the
problem with boundary conditions as a linear combination of them, solving for
the coefficients that minimize the error on the boundary. Typically, the
choices have been rational functions [28] or products of Bessel functions and
trigonometric polynomials centered at certain points. The different choices of
these functions have a big impact on the performance of the method. Recently,
Gopal and Trefethen [24, 25] have developed a new way of selecting the base
functions in such a way as to yield root exponential convergence (the
lightning Laplace solver). We stress that these methods produce accurate
approximations but there is no explicit control of the error with respect to
the true solution. This is handled a posteriori with a perturbative analysis
of the approximations.

The use of computers to prove mathematically rigorous theorems has become
increasingly popular in the last 20 years and many goals and theories have
been developed in this blooming field. Floating point arithmetic errors are
handled and controlled via interval arithmetic, where real numbers are
replaced by real intervals and all the errors are propagated throughout the
calculations. We refer to the book [53] for an introduction to validated
numerics, and to the survey [22] and the recent book [43] for a more specific
treatment of computer-assisted proofs in PDE. We also mention the work of
Tanaka [52] where he also controls the nodal line of an elliptic problem (as
opposed to an eigenvalue problem in our case) using computer-assisted
techniques, and Dahne–Salvy [13] enclosing eigenvalues of the Laplacian on
spherical triangles with techniques very similar to those used here.

The paper is organized as follows. In Section 2 we present the candidate for
the counterexample and discuss how it was found. In Section 3 we explain the
separation of the first four eigenvalues from the rest. Section 4 gives
details of how the approximate eigenfunctions are constructed and in Section 5
we make use of these approximations for isolating the second eigenvalue.
Finally in Section 6 we prove that the nodal line of the second eigenfunction
is closed and conclude the proof of Theorem 1.2. Section 7 contains details
about the implementation of the computer assisted parts.

##  2 Finding a candidate

(a) The candidate domain is a hexagon with side length $1$ with holes that are
equilateral triangles of height $\frac{6}{27}$ with centers placed at a
distance of $\frac{11}{27}$ from the origin.

(b) Plot highlighting the nodal line of $\tilde{u}_{2}$, generated by plotting
$\operatorname{sign}(\tilde{u}_{2})\log(|\tilde{u}_{2}|)$.

Figure 1: The counterexample and an approximation of its nodal domains.

The main idea behind the choice of the domain was to start from Hoffmann-
Ostenhof–Hoffman-Ostenhof–Nadirashvili’s construction from a disk and carve as
few (but possibly large) holes as possible. We tried to work with domains as
symmetric as possible and holes with few sides to reduce the computational
cost. However, due to the lower bound of Theorem 3.1 being restricted to
polygonal domains we chose the domain to be of polygonal shape. We found many
instances of domains for which the nodal line was closed, though the problem
seems to be quite sensitive to the position and shape of the holes, and small
perturbations destroy the closedness of the nodal line due to the very small
nature of the relevant numbers, see Figure 2 for a few different candidates.
In the end we settled for a domain given by a hexagon with six holes in it.
The side length of the hexagon is normalized to 1 and the holes are
equilateral triangles of height $\frac{6}{27}$ with their centers placed a
distance $\frac{14}{27}$ from the origin, see Figure 1(a). The search was done
using the PDE Toolbox in Matlab and the first five approximate eigenvalues for
the final domain were

$\lambda_{1}=31.0432,\ \lambda_{2}=63.2104,\ \lambda_{3}=63.7259,\
\lambda_{4}=63.7259,\ \lambda_{5}=68.2629.$ (2)

Figure 2: Four domains for which the nodal line seems to be closed. Similarly
to Figure 1(b) we plot $\operatorname{sign}(u)\log(|u|)$. The rightmost one is
the one we finally chose and the only one for which we have proved that the
nodal line indeed is closed.

##  3 Separating the first four eigenvalues

The first four eigenvalues will be separated from the rest by computing a
lower bound for the fifth eigenvalue. We make use of recent results by [38].
The procedure is exactly the same as in [23].

The starting point is a triangulation of the domain as given in Figure 3(a).
It consists mostly of equilateral triangles except next to the boundary of the
holes where the triangles are cut in half (these are colored red in the
figure). The basis functions are indexed by the interior edges of the
triangulation: For an edge $E$ between two triangles $\tau_{1}$ and $\tau_{2}$
the basis function $\psi_{E}$ is the unique function supported on
$\tau_{1}\cup\tau_{2}$ such that the restriction to each triangle is affine,
takes the value 1 at the midpoint of $E$ and the value 0 at the midpoints of
the other edges of $\tau_{1}$ and $\tau_{2}$. In our case there are three
kinds of edges, type 1 where $\tau_{1}$ and $\tau_{2}$ are both equilateral
(these correspond to the majority of edges), type 2 where only one of them is
equilateral and type 3 where none of them are equilateral.

(a) Triangulation of the domain.

(b) Triangulation of the domain around the holes.

Figure 3: Meshing of the counterexample.

The weak formulation of the problem (1) reads

$\int_{\Omega}\nabla u\cdot\nabla\psi=\lambda\int_{\Omega}u\psi.$

We define the coefficients of the stiffness and mass matrices
$A=(a_{EF}),B=(b_{EF})$ by

$a_{EF}=\int_{\Omega}\nabla\psi_{E}\cdot\nabla\psi_{F},\quad
b_{EF}=\int_{\Omega}\psi_{E}\psi_{F}$

which leave us with solving discrete system

$Ax=\lambda Bx.$ (3)

Here the vector $x=(x_{E})$ corresponds to the discrete solution $u=\sum
x_{E}\psi_{E}$. For our choice of triangulation the mass matrix will be
diagonal with coefficients

$\frac{h^{2}}{2\sqrt{3}},\quad\frac{3h^{2}}{8\sqrt{3}},\quad\frac{h^{2}}{4\sqrt{3}},\quad$

for edges of type 1, 2 and 3 respectively. Here $h$ is the side length of the
equilateral triangles. The matrix $A$ is sparse and on the diagonal the
coefficients are

$\frac{8}{\sqrt{3}},\quad\frac{12}{\sqrt{3}},\quad\frac{4}{\sqrt{3}},\quad$

for edges of type 1, 2 and 3 respectively. The off-diagonal entries are zero
if the two edges do not have a common triangle or $\frac{-2}{\sqrt{3}}$ if
they do. Since $B$ is diagonal it is easily inverted, which allows us to
reduce the problem to solving the matrix eigenvalue problem

$Mx=\lambda x,\quad M=B^{-1}A.$ (4)

The solution of the finite element problem can then be linked to the
continuous problem using the following Theorem:

######  Theorem 3.1.

[38, Theorem 2.1, Remark 2.2] Consider a polygonal domain $\Omega$ with a
triangulation so that each triangle has diameter at most $h$. Let
$\lambda_{k}$ be the $k$-th solution of the eigenvalue problem (1) in $\Omega$
and $\lambda_{h,k}$ the k-th eigenvalue of the Crouzeix-Raviart discretized
problem (4) in $\Omega$. Then

$\frac{\lambda_{h,k}}{1+C_{h}^{2}\lambda_{h,k}}\leq\lambda_{k},$ (5)

where $C_{h}\leq 0.1893h$ is a constant.

Due to the monotonicity of (5) this reduces the problem of lower bounding the
fifth eigenvalue of the domain: $\lambda_{5}$, to lower bounding the fifth
eigenvalue of (4): $\lambda_{h,5}$.

To get a lower bound for $\lambda_{h,5}$ we will follow the same procedure as
in [23], making use of Gershgorin’s generalized disks theorem [19]. Before we
can do that we need the following lemma:

######  Lemma 3.2.

[23, Lemma 2.4] Let $v_{1},\dots,v_{m}$ be vectors in $\mathbb{R}^{m}$ and
$s>0$ such that $|\langle v_{i},v_{j}\rangle-\delta_{ij}|\leq s$ and suppose
that $8ms<1$. Then there exists an orthonormal set of vectors
$w_{1},\dots,w_{m}\in\mathbb{R}^{m}$ such that $\|v_{i}-w_{i}\|\leq\sqrt{3s}$.

We will use it in the following way: we first find an approximate orthonormal
basis of eigenvectors $\\{v_{i}\\}$ of $M$, this is easily done using standard
eigenvalue routines. Let $\tilde{Q}$ be the matrix that has them as columns,
and compute an enclosure for the almost-diagonal matrix
$\tilde{D}=\tilde{Q}^{T}M\tilde{Q}$. Using the lemma we can find an
orthonormal matrix $Q$ which is close to $\tilde{Q}$. Let $D=Q^{T}MQ$, which
is not necessarily diagonal but has the same eigenvalues as $M$ due to the
orthogonality of $Q$. We can obtain rigorous enclosures of the entries of $D$
using $\tilde{D}$ in the following way

$\displaystyle|D_{ij}-\tilde{D}_{ij}|$ $\displaystyle=|\langle
w_{i},Mw_{j}\rangle-\langle v_{i},Mv_{j}\rangle|$ $\displaystyle\leq|\langle
w_{j}-v_{j},Mv_{i}\rangle+\langle w_{j}-v_{j},M(w_{i}-v_{i})\rangle+\langle
w_{i}-v_{i},Mv_{j}\rangle|$
$\displaystyle\leq\sqrt{3s}(\|Mv_{i}\|+\|Mv_{j}\|)+4s\|M\|_{2},$

with $s$ is as in the lemma and using the symmetry of $M$. Observe that
$\|Mv_{i}\|$ can be computed explicitly and the upper bound
$\|M\|_{2}\leq\|M\|_{\text{Frob}}$ is easily computed.

Finally applying Gershgorin’s generalized disks theorem to the matrix $D$, of
which we have sharp bounds, we can separate the spectrum of $M$ in two
components, one of which will contain the first 4 eigenvalues and the other of
which will contain the rest. The lower bound we get for the second component
is $66.2862$ which is thus also a lower bound for $\lambda_{h,5}$. A direct
application of Theorem 3.1 then gives us the lower bound $66.0709$ for
$\lambda_{5}$. Note that this is well above the approximate value for
$\lambda_{4}$ computed in Section 2.

##  4 Constructing approximate solutions

In this section we explain how to compute approximations for the first four
eigenfunctions which we will then use to isolate the second eigenfunction and
analyse its nodal line. We will make use of the Method of Particular Solutions
(MPS) to the compute the approximations. The version of MPS that we give here
is due to Betcke and Trefethen [9, 8]. One starts writing the eigenfunction as
a linear combination of functions $\phi_{i}\,(1\leq i\leq N)$ that satisfy the
equation $-\Delta\phi_{i}=\lambda\phi_{i}$ in the domain but with no boundary
conditions, where $\lambda$ is taken as a parameter and is part of the
problem. The coefficients are then chosen to approximate the boundary
condition we want $u$ to satisfy. In the case of a zero Dirichlet boundary
condition this amounts to finding $\lambda$ and a non-zero linear combination
for which the boundary values are as close to zero as possible. The linear
combination is determined by fixing $\lambda$ and taking $m_{b}$ collocation
points on the boundary, $\\{x_{k}\\}_{k=1}^{m_{b}}\subset\partial\Omega$, then
choosing the linear combination to minimize its values on the collocation
points in the least squares sense. This alone is not quite enough, increasing
the number $N$ of elements in the basis leads to existence of linear
combinations very close to 0 inside the domain $\Omega$. The version by Betcke
and Trefethen handles this by also adding a number $m_{i}$ of interior points,
$\\{y_{l}\\}_{l=1}^{m_{i}}\subset\Omega$, and taking the linear combination to
stay close to unit norm on these. This is accomplished by considering the two
matrices $A_{B}=(\phi_{i}(x_{k}))$ and $A_{I}=(\phi_{i}(y_{l}))$ which are
combined into a matrix whose _QR_ factorization gives an orthonormal basis of
these function evaluations

$A=\begin{bmatrix}A_{B}\\\ A_{I}\end{bmatrix}=\begin{bmatrix}Q_{B}\\\
Q_{I}\end{bmatrix}R=:QR.$

The right singular vector $v$ corresponding to the smallest singular value
$\sigma=\sigma(\lambda)$ of $Q_{B}$ for a given $\lambda$ is a good candidate
for the eigenfunction when $\sigma(\lambda)$ is small. In our case we do not
require very high precision for the first, third and fourth eigenvalue and the
approximate values from (2) gives good enough approximate eigenfunctions. The
second eigenfunction needs slightly higher precision and this approximate
value is not good enough. Instead we search for a $\lambda$ around it which
minimizes $\sigma(\lambda)$ using Brent’s method. This minimizer is the value
given in Table 1.

We will use three types of basis functions $\phi_{i}$, all of them products of
Bessel functions and trigonometric functions. They are all given in polar
coordinates centered around a certain point (which may be different from basis
function to basis function). The first type is the one used in the original
version of MPS, it is centered around a vertex of the domain. If the angle of
the vertex is $\pi/\alpha$ they take the form

$\phi_{\alpha,k}=J_{\alpha k}(\sqrt{\lambda}r)\sin\alpha k\theta$ (6)

in polar coordinates centered around the vertex and $\theta=0$ taken along one
of the boundary segments. Unless $\alpha$ is an integer these functions have a
branch cut in $\theta$ and are therefore only suitable for vertices where the
branch cut can be placed outside of the domain, in particular they are not
suited for placement at the vertices of the holes of our domain. The two other
types of functions were recently introduced by Gopal and Trefethen [25]. One
is centered around points (called charges) that accumulate root-exponentially
near a vertex of the domain and take the form

$\phi_{\text{ext}}(r,\theta)=Y_{0}\left(r\sqrt{\lambda}\right),\qquad\phi_{\text{ext}}^{\text{c}}(r,\theta)=Y_{1}\left(r\sqrt{\lambda}\right)\cos\theta,\qquad\phi_{\text{ext}}^{\text{s}}(r,\theta)=Y_{1}\left(r\sqrt{\lambda}\right)\sin\theta.$
(7)

In this case $\theta=0$ is taken along the bisector of the vertex. The final
type is an interior expansion which in our case is placed at the center of the
domain with $\theta=0$ taken along the positive $x$-axis

$\phi_{0}(r,\theta)=J_{0}\left(r\sqrt{\lambda}\right),\qquad\phi_{j}^{\text{c}}(r,\theta)=J_{j}\left(r\sqrt{\lambda}\right)\cos
j\theta,\qquad\phi_{j}^{\text{s}}(r,\theta)=J_{j}\left(r\sqrt{\lambda}\right)\sin
j\theta.$ (8)

Here $J$ and $Y$ are the Bessel functions of the first and second kind
respectively.

The basis functions will be used in the following way:

  1. 1.

At each vertex of the hexagon we place $2N_{1}$ basis functions of the first
type (6), $\phi_{\alpha,k}$ ($1\leq k\leq 2N_{1}$). The interior angle in this
case is $2\pi/3$ so $\alpha=3/2$. This gives $2N_{1}$ free coefficients for
each vertex for a total of $12N_{1}$ free coefficients.

  2. 2.

At each vertex of the holes we put $N_{2}$ charges and basis functions of the
second kind (7) around those. This gives $3N_{2}$ free coefficients for each
vertex for a total of $54N_{2}$ charges counting all the vertices.

  3. 3.

In the center of the domain we place basis functions of the third type (8)
with $0\leq j\leq 6N_{3}-1$, for a total of $12N_{3}-1$ free coefficients
($j=0$ giving only one free coefficient).

In the computations we fix $n$ and let $N_{1}=N_{3}=n$ and $N_{2}=3n$. The
number of free coefficients is then $12n+54\cdot 3n+12n-1=186n-1$. See Figure
4 for a schematic of the placement of the basis functions and charges.

Figure 4: Schematic of placement of charges. The small dots indicate placement
of charges for the basis functions of the second type, only a few charges are
shown here. For the first and second eigenfunction the expansions with the
same color can be taken to have the same coefficients due to symmetry. For the
third and fourth eigenfunction the symmetries are different.

We can reduce the number of free coefficients substantially by making use of
the symmetries of the domain and the eigenfunctions. We expect the first and
second eigenfunction to have 6-fold symmetry and also be even with respect to
the $x$-axis. This allows us to reduce the number of free coefficients by both
fixing some of the expansions to have the same coefficients and by avoiding
the use of some terms in the expansions which do not satisfy the required
symmetry.

  1. 1.

For the expansions at the vertices of the hexagon we require that they all
have the same coefficients, reducing the number of free coefficients to
$2N_{1}$ instead of $12N_{1}$. Furthermore, due to their evenness it is enough
to consider only even values of $k$, further reducing the number of free
coefficients to $N_{1}$.

  2. 2.

For the expansions at the vertices of the holes the 6-fold symmetry allows us
to reduce the number of free coefficients to $9N_{2}$. We can reduce it
further by using that it has to be even. For the expansion at the outer vertex
of each hole we only have to consider $\phi_{\text{ext}}(r,\theta)$ and
$\phi_{\text{ext}}^{\text{c}}(r,\theta)$, so the number of coefficients for
that vertex is reduced to $2N_{2}$. Finally we can notice that the expansions
at the two inner vertices of the holes are symmetric due to the evenness, so
must also have the same coefficients. This reduces the total number of free
coefficients to $5N_{2}$.

  3. 3.

The inner expansion should be even and have 6-fold symmetry. This means we can
skip the use of $\phi_{j}^{\text{s}}(r,\theta)$ and only have to consider
$j=6i$ with $0\leq i\leq N_{3}-1$. Reducing the number of free coefficients to
$N_{3}$.

This reduces the number of free coefficients from $184n$ to $17n$, more than a
10x improvement.

We expect that the third eigenvalue is even with respect to the $x$-axis and
odd with respect to the $y$-axis, whereas the fourth one will satisfy the
opposite symmetries. Using these symmetries we take $45n$ free coefficients in
the first case and $42n$ in the second case.

The amount of accuracy we need is different for the different eigenfunctions.
For the first eigenfunction we only need enough accuracy to separate
$\lambda_{1}$ from $\lambda_{2}$, and since they are far away from each other
very low accuracy is enough. For the third and fourth eigenfunction we need
slightly higher accuracy since $\lambda_{3}$ and $\lambda_{4}$ are much closer
to $\lambda_{2}$, we also need slightly higher accuracy to be able to handle
the fact that $\lambda_{3}$ and $\lambda_{4}$ presumably correspond to a
double eigenvalue. The second eigenfunction is the most challenging one: we
need much higher accuracy to be able to isolate the nodal line. Details of the
computations are given in Table 1.

Eigenfunction | Eigenvalue | $n$ | Free coefficients | Collocation points
---|---|---|---|---
$\tilde{u}_{1}$ | $\tilde{\lambda}_{1}=31.0432$ | 1 | 17 | 51
$\tilde{u}_{2}$ | $\tilde{\lambda}_{2}=63.20833598626884$ | 28 | 476 | 7616
$\tilde{u}_{3}$ | $\tilde{\lambda}_{3}=63.7259$ | 6 | 270 | 2160
$\tilde{u}_{4}$ | $\tilde{\lambda}_{4}=63.7259$ | 6 | 252 | 2016
Table 1: Details about computation of the approximate eigenpairs.

We finish this section with four approximate eigenpairs
$(\tilde{\lambda}_{1},\tilde{u}_{1})$, $(\tilde{\lambda}_{2},\tilde{u}_{2})$,
$(\tilde{\lambda}_{3},\tilde{u}_{3})$ and
$(\tilde{\lambda}_{4},\tilde{u}_{4})$, as detailed in Table 1 and Figure 5. In
the next two sections we will use these approximations to isolate the second
eigenvalue and prove that the nodal line of the second eigenfunction is
closed.

Figure 5: The four approximate eigenfunctions, starting with $\tilde{u}_{1}$
in the top left and ending with $\tilde{u}_{4}$ in the bottom right.

##  5 Isolating the second eigenvalue

The main tool in this section is the following theorem by Fox, Henrici and
Moler from the original MPS article

######  Theorem 5.1.

[17, 42] Let $\Omega\subset\mathbb{R}^{n}$ be bounded. Let $\tilde{\lambda}$
and $\tilde{u}$ be an approximate eigenvalue and eigenfunction—that is, they
satisfy $\Delta\tilde{u}+\tilde{\lambda}\tilde{u}=0$ in $\Omega$ but not
necessarily $\tilde{u}=0$ on $\partial\Omega$. Define

$\mu=\frac{\sqrt{|\Omega|}\sup_{x\in\partial\Omega}|\tilde{u}(x)|}{\|\tilde{u}\|_{2}}.$
(9)

where $|\Omega|$ is the area of the domain. Then there exists an eigenvalue
$\lambda$ such that

$\frac{|\tilde{\lambda}-\lambda|}{\lambda}\leq\mu.$ (10)

This allows us to calculate enclosures of the approximate eigenvalues computed
in the previous section. To do that we need to find an upper bound for $\mu$
in (9), which in turn means we need an upper bound for the maximum on the
boundary and a lower bound for the $L^{2}$ norm of the approximate
eigenfunctions.

The upper bound for the maximum is computed using a combination of interval
arithmetic and Taylor expansions. The sides are divided into small segments
and on the midpoint of each segment a Taylor polynomial of the approximate
eigenfunction is computed. The maximum of the polynomial is then bounded using
classical interval arithmetic and the error term for the polynomial is added.
For more details see [13, Section 2.1]. Thanks to the symmetries satisfied by
the approximate eigenfunctions we only have to check the maximum on parts of
the boundary.

For lower bounding the norm we use the same procedure as in [23] and [13].
Consider a subset of the domain, $\Omega^{\prime}\subset\Omega$. If $u$ does
not vanish in $\Omega^{\prime}$ then without loss of generality it can be
assumed to be positive there and then, since $-\Delta u=\lambda u>0$, $u$ is
superharmonic in $\Omega^{\prime}$ and satisfies
$\inf_{\Omega^{\prime}}u\geq\inf_{\partial\Omega^{\prime}}u$. Thus a lower
bound for $|u|$ on $\partial\Omega^{\prime}$ yields a lower bound for $|u|$
inside $\Omega^{\prime}$. To determine that $u$ does not vanish on
$\Omega^{\prime}$ a lower bound for $|u|$ on $\partial\Omega^{\prime}$ is
computed with the same techniques as when upper bounding the maximum. Once it
is determined that $u$ has a fixed sign on $\partial\Omega^{\prime}$ (which we
assume to be positive) the key observation is that $u$ cannot be negative
inside $\Omega^{\prime}$ if $\Omega^{\prime}$ is small enough. Indeed, if
$\Omega^{\prime\prime}\subset\Omega^{\prime}$ is a maximal domain where $u<0$,
then $u=0$ on $\partial\Omega^{\prime\prime}$ and thus $\lambda$ is an
eigenvalue for $\Omega^{\prime\prime}$. If the area of $\Omega^{\prime}$ is
small enough this scenario can be ruled out using the Faber-Krahn inequality.
Details about the choice of $\Omega^{\prime}$ for the different eigenfunctions
are given in Section 7.

Upper bounds for $\mu$ together with upper bounds for the maximum on the
boundary and enclosures of the eigenvalues are given in Table 2. From Section
3 we know that there are at most four eigenvalues below $66.0709$. If all of
the four enclosures were isolated we could have concluded that we had isolated
all four eigenvalues and we would know their indices in the spectrum. However
the enclosures coming from $\tilde{u}_{3}$ and $\tilde{u}_{4}$ overlap and we
have to handle that.

Eigenfunction | $\mu\leq$ | $\sup_{x\in\partial\Omega}\leq$ | Enclosure for eigenvalue
---|---|---|---
$\tilde{u}_{1}$ | 0.14 | 0.012 | $[30\pm 6.1]$
$\tilde{u}_{2}$ | $8.26\cdot 10^{-5}$ | $1.01\cdot 10^{-6}$ | $[63.21\pm 6.89\cdot 10^{-3}]$
$\tilde{u}_{3}$ | 0.00186 | $2.68\cdot 10^{-5}$ | $[64\pm 0.393]$
$\tilde{u}_{4}$ | 0.00215 | $3.1\cdot 10^{-5}$ | $[64\pm 0.411]$
Table 2: Upper bounds of $\mu$ and the value on the boundary of the four
approximate eigenfunctions together with computed enclosures for the
corresponding eigenvalue.

###  5.1 Handling the 2-cluster

Since the enclosures coming from $\tilde{\lambda}_{3}$ and
$\tilde{\lambda}_{4}$ overlap we cannot be sure that they indeed correspond to
two different eigenvalues. From the plots of the corresponding eigenfunctions
in Figure 5 it does indeed seem like they do correspond to different
eigenvalues but what we will prove is a slightly weaker statement which is
enough for what we want to do. We will prove that there are at least two
eigenvalues in an interval slightly larger than the two enclosures.

The proof will be based on the fact that $\tilde{u}_{3}$ and $\tilde{u}_{4}$
are not proportional to each other so as long as the error bounds for them are
sufficiently small this would imply that the corresponding exact solutions can
not possibly correspond to the same eigenfunction. The error bounds we will
use are given by the following theorem:

######  Theorem 5.2.

[42, Theorem 3] Let $\Omega\subset\mathbb{R}^{2}$ be bounded. Let
$\tilde{\lambda}$ and $\tilde{u}$ be an approximate eigenvalue and
eigenfunction—that is, they satisfy
$\Delta\tilde{u}+\tilde{\lambda}\tilde{u}=0$ in $\Omega$ but not necessarily
$\tilde{u}=0$ on $\partial\Omega$. Let

$\mu=\frac{\sqrt{|\Omega|}\sup_{x\in\partial\Omega}|\tilde{u}(x)|}{\|\tilde{u}\|_{2}}.$

Then there is an eigenvalue $\lambda_{k}$ satisfying the same bounds as in
Theorem 5.1 and a corresponding eigenfunction $u_{k}$. Let

$g(x)=\left(\int_{\Omega}G(x,y)^{2}dy\right)^{1/2}$

where $G(x,y)$ is the Green’s function and

$\alpha=\min_{\lambda_{n}\not=\lambda_{k}}\frac{|\lambda_{n}-\tilde{\lambda}|}{|\lambda_{n}|}.$

Then for any $x\in\Omega$ we have the following bound for $\tilde{u}$

$|\tilde{u}(x)-u_{k}(x)|\leq\left(\sup_{x\in\partial\Omega}|\tilde{u}(x)|\right)\left(1+g(x)\tilde{\lambda}\left(\frac{1}{1-\mu}+\frac{1}{\alpha}\left(1+\frac{\mu^{2}}{\alpha^{2}}\right)\right)\right).$

To compute this bound we need an upper bound for $g(x)$ and a lower bound for
$\alpha$. Using [3, Theorem 2.4, Corollary 2.2] with $p=2$, we obtain that

$\displaystyle g(x)\leq\frac{1}{4\pi}\sqrt{2|\Omega|}$ (11)

######  Remark 5.3.

The bound (11) is far from optimal (especially when $x\to\partial\Omega$), but
we have preferred to keep a simple, uniform bound.

######  Remark 5.4.

It may be possible to obtain $L^{\infty}$ bounds by deriving higher order
estimates and using the Sobolev embedding as in [47] but our method does not
require any additional explicit constant.

We are now ready to handle the 2-cluster. Let $\Lambda$ be the union of the
enclosures of the eigenvalues $\tilde{\lambda}_{3}$ and $\tilde{\lambda}_{4}$
as given by Theorem 5.1. Then we have the following result.

######  Lemma 5.5.

Let $r$ be the radius of $\Lambda$ and let $\Lambda^{\prime}$ be the interval
with the same midpoint as $\Lambda$ and radius $\frac{17r}{16}$. Then there
are at least two eigenvalues in the interval $\Lambda^{\prime}$.

###### Proof.

The proof is by contradiction. Assume that there is only one eigenvalue in
$\Lambda^{\prime}$. From Theorem 5.2 we get that there are two eigenfunctions
$u_{3}$ and $u_{4}$ with eigenvalues $\lambda_{3}$ and $\lambda_{4}$ in
$\Lambda$. By our assumption, we must have $\lambda_{3}=\lambda_{4}$ and also
$u_{3}=Cu_{4}$ for some $C\in\mathbb{R}$.

Since $\lambda_{3}=\lambda_{4}$ lie in $\Lambda$ and we assume that there are
no other eigenvalues in $\Lambda^{\prime}$ we get that $\alpha$ in Theorem 5.2
is lower bounded by $\frac{r}{16}$. With a lower bound for $\alpha$ and upper
bounds for $g(x)$, $\sup_{x\in\partial\Omega}|\tilde{u}(x)|$ and $\mu$ we can
use Theorem 5.2 to compute an upper bound $d_{3}$ for
$|\tilde{u}_{3}(x)-u_{3}(x)|$ and $d_{4}$ for $|\tilde{u}_{4}(x)-u_{4}(x)|$.

Now consider the points $p_{1}=(1/2,1/2)$ and $p_{2}=(-1/2,1/2)$. We can
evaluate $\tilde{u}_{3}$ and $\tilde{u}_{4}$ at these points. We find that
$\tilde{u}_{3}(p_{1})+d_{3}<0\text{ and }\tilde{u}_{4}(p_{1})+d_{4}<0$ so both
$u_{3}$ and $u_{4}$ must have the same sign at $p_{1}$, this means that
$u_{3}=Cu_{4}$ for some $C>0$ and in particular $u_{3}$ and $u_{4}$ must have
the same sign everywhere. Furthermore we find that
$\tilde{u}_{3}(p_{2})-d_{3}>0$ so $u_{3}$ is positive at $p_{2}$ but
$\tilde{u}_{4}(p_{2})+d_{3}<0$ so $u_{4}$ must be negative at $p_{2}$. This
would contradict that $u_{3}$ and $u_{4}$ have the same sign and hence there
must be at least two eigenvalues in $\Lambda^{\prime}$. ∎

######  Remark 5.6.

While Lemma 5.5 proves that there are at least two eigenvalues in
$\Lambda^{\prime}$ it does not prove that there has to be a double eigenvalue.

######  Remark 5.7.

We remark that it is fundamental for the cluster to be a 2-cluster in the
proof. Otherwise, there are no lower bounds on $\alpha$ available, in order to
apply Theorem 5.2. We outline a different (though more costly) strategy for
the general case.

The idea is to make use of the symmetries of the domain and the approximate
eigenfunctions and compute bounds on the smaller domain
$\Omega^{\prime}\subset\Omega$ which are then extended to $\Omega$. For
example, for $\tilde{u}_{3}$ we could do the following. Consider
$\tilde{u}_{3}$ as an approximate eigenfunction on the domain
$\Omega^{\prime}$ given by the right half of $\Omega$. If we can compute
$L^{\infty}$ bounds for $\tilde{u}_{3}$ on $\Omega^{\prime}$ then these bounds
also apply to $\Omega$ extending by symmetry. When computing bounds on
$\Omega^{\prime}$ using Theorem 5.2 we would get exactly the same bound for
$\mu$ since by construction, $\tilde{u}_{3}$ is identically equal to zero on
the only additional boundary and the factor difference in norm and area scale
out. In $\Omega^{\prime}$, $\tilde{u}_{3}$ corresponds to a simple eigenvalue
and we could get a lower bound for $\alpha$ without having to deal with a
cluster of eigenvalues. This would require more work than Lemma 5.5 since we
would have to control the spectrum of $\Omega^{\prime}$ as well, but it has
the benefit that it would work for clusters with more than two eigenvalues.

We have now proved that there must be one eigenvalue in $[30\pm 6.1]$, one in
$[63.21\pm 6.89\cdot 10^{-3}]$ and two in $\Lambda^{\prime}$. Since we have
already proved that there are at most four eigenvalues below $66.0709$ we can
conclude that we have found them all. In particular this means that
$(\tilde{\lambda}_{2},\tilde{u}_{2})$ must indeed correspond to the second
eigenpair and the eigenvalue closest to it lies in $\Lambda^{\prime}$.

##  6 Analysing the nodal line and conclusion of the proof

In this section we will assume that the sign of $\tilde{u}_{2}$ is taken such
that it is positive at the center. To prove that the nodal line is isolated we
will construct a closed path $\Gamma$ around the center of the domain. Let
$\tilde{\Omega}$ be the region enclosed by $\Gamma$, $\Gamma$ will be taken
such that $\tilde{\Omega}$ does not intersect the boundary of $\Omega$. We
will prove that $u_{2}$ is strictly negative on $\Gamma$ and positive in some
point in $\tilde{\Omega}$. This would imply that $u_{2}$ changes sign
somewhere inside of $\tilde{\Omega}$ and hence at least part of the nodal line
must be inside it. Since $u_{2}$ is strictly negative on $\Gamma$ we also get
that the nodal line can not cross $\Gamma$, thus there is a closed curve in
$\tilde{\Omega}$ belonging to the nodal line. By Courant’s nodal domain
Theorem [12] there can only be two connected components and hence the whole
nodal line is closed and contained in $\tilde{\Omega}$, therefore not touching
$\partial\Omega$.

The choice of $\Gamma$ is seen in Figure 6. It consists of a straight line
between $(d,0)$ and $(d,\tan(\pi/6)d)$ (the solid part of the red line in
Figure 6) which is mirrored in the $x$-axis and extended in a 6-fold symmetry
way. Since $\tilde{u}_{2}$ satisfies the same symmetries it is enough to bound
it on this straight line to get a bound on $\Gamma$. The value of $d$ is
chosen so that the value of $\tilde{u}_{2}(d,0)$ is as negative as possible.
We picked $d=0.38483177115481165$.

Figure 6: Shows the path $\Gamma$, enclosing the region $\tilde{\Omega}$, in
red. By symmetry arguments it is enough to consider the value of
$\tilde{u}_{2}$ on the solid part of the line.

We compute an upper bound for $\tilde{u}_{2}$ on $\Gamma$ using the same
techniques as when computing the bound on the boundary. In this case it is
much less costly since there are less cancellations between the terms of
$\tilde{u}_{2}$ and we are further away from the charge points. We get the
upper bound $\gamma=-4.4929\cdot 10^{-5}$. We are left to prove that the error
bound for $\tilde{u}_{2}$ from Theorem 5.2 is less than $|\gamma|$ to conclude
that $u_{2}$ must be negative on all of $\Gamma$.

From Lemma 5.5 we know that the closest eigenvalue to $\tilde{\lambda}_{2}$ is
in the interval $\Lambda^{\prime}$. This gives us the bound $\alpha\geq
0.3727$. Together with the upper bound for $g(x)$ from (11) and the upper
bounds for $\mu$ and $\max_{x\in\partial\Omega}|\tilde{u}_{2}(x)|$ from Table
2, Theorem 5.2 gives us

$|u_{2}(x)-\tilde{u}_{2}(x)|\leq 4.2162\cdot 10^{-5}.$

Since $|\gamma|>4.2162\cdot 10^{-5}$ we conclude that $u_{2}$ is negative on
all of $\Gamma$.

Finally we have to prove that $u_{2}$ is positive at some point in
$\tilde{\Omega}$, for that we just take the point $(1/10,0)$ for which we have
$\tilde{u}_{2}(1/10,0)\in[0.01342\pm 3.27\cdot 10^{-6}]$. Since this is
greater than the error bound for $\tilde{u}_{2}$ we can conclude that $u_{2}$
is positive on at least one point in $\tilde{\Omega}$, which implies that the
nodal line for $u_{2}$ must be fully contained in $\tilde{\Omega}$ and hence
is closed. This finishes the proof of Theorem 1.2.

##  7 Details of the Implementation

The code111Available at https://github.com/Joel-Dahne/PaynePolygon.jl. for the
computer assisted parts is implemented in Julia [10] and makes use of Arb [32]
for the rigorous parts of the numerics. The part related to Section 3 is
implemented from scratch whereas the rest relies heavily on the code from [13]
with added support for new geometries and new types of basis functions.
Roughly the code is divided into three parts:

  1. 1.

The code for computing the lower bound of $\lambda_{h,5}$ using the FEM
method. The computation of $Q$ uses standard eigenvalue routines and the
verification is then easily done using Arb. It is complicated by the fact that
the matrix $M$ is very large ($6084\times 6084$) and that standard double
precision computations give an $s$ in Lemma 3.2 which is not quite small
enough for the separation to work. To get higher precision we do the
computations with so called double-doubles [50] through David K. Zhang’s Julia
library for multi-float
computations222https://github.com/dzhang314/MultiFloats.jl. This has the
drawback that there are less specialised eigenvalue routines and the
computations thus take longer, hours instead of minutes, but otherwise it
works well.

  2. 2.

The code for computing the approximate eigenfunctions is mostly the same as in
[13] with some modifications for the planar case. The collocation points are
placed root-exponentially close to the vertices, in line with the
recommendations of [25]. For the lightning charges we use standard double
precision. Between the expansions at the vertices of the hexagon and the
expansion at the center we see very large cancellations and we found that
increasing the precision helped. For the first eigenfunction 64 bits was
enough and for the third and fourth we found 128 bits to suffice. For the
third eigenfunction the cancellations were even larger and we ended up using
384 bits of precision.

  3. 3.

The computations required for the verification of the approximate solutions,
mostly the lower bound of the norm and the upper bound on the maximum for the
approximate eigenfunctions, is for the most part exactly the same as in [13].
The subset of the domain where the norm is lower bounded is tuned for each
eigenfunction and given in Figure 7.

The calculations were run on relatively old hardware, the FEM part on an Intel
Xeon CPU E5-2620 with 48GB of memory and the rest on an Intel Core i7-3770
with 16GB of memory. The most time consuming part of the computations is the
FEM method. The computation of $\tilde{Q}$ took around 30 hours and then an
additional 2 hours for the verification with Arb. Computation of the
approximate eigenfunctions took around 10 seconds for $\tilde{u}_{1}$, 6 hours
for $\tilde{u}_{2}$ and 5 minutes each for $\tilde{u}_{3}$ and
$\tilde{u}_{4}$. The bounds for the norm took less than 20 minutes in total
whereas the bound for the maximum took around 10 seconds for $\tilde{u}_{1}$,
6 hours for $\tilde{u}_{2}$ and half an hour each for $\tilde{u}_{3}$ and
$\tilde{u}_{4}$.

Figure 7: The part of the domain where the norm of each eigenfunction is lower
bounded is marked with red lines. For symmetry reasons it is enough to bound
the eigenfunctions on the parts with solid lines.

## Acknowledgments

JD was partially supported by the European Research Council through ERC-
StG-852741-CAPA. JGS was partially supported by the European Research Council
through ERC-StG-852741-CAPA, by NSF through Grant NSF DMS-1763356 and by the
Princeton Summer Program for Mathematics Majors. KH was partially supported by
the Princeton Summer Program for Mathematics Majors. We thank Uppsala
University for computing facilities (Haddock Cluster).

## References

  * [1] Alessandrini, G. Nodal lines of eigenfunctions of the fixed membrane problem in general convex domains.  Comment. Math. Helv. 69, 1 (1994), 142–154.
  * [2] Antunes, P. R. S., and Valtchev, S. S. A meshfree numerical method for acoustic wave propagation problems in planar domains with corners and cracks.  J. Comput. Appl. Math. 234, 9 (2010), 2646–2662.
  * [3] Bandle, C. Isoperimetric inequalities and applications, vol. 7 of Monographs and Studies in Mathematics.  Pitman (Advanced Publishing Program), Boston, Mass.-London, 1980.
  * [4] Barnett, A. H., and Hassell, A. Boundary quasi-orthogonality and sharp inclusion bounds for large Dirichlet eigenvalues.  SIAM J. Numer. Anal. 49, 3 (2011), 1046–1063.
  * [5] Beattie, C., and Goerisch, F. Methods for computing lower bounds to eigenvalues of self-adjoint operators.  Numer. Math. 72, 2 (1995), 143–172.
  * [6] Beck, T., Canzani, Y., and Marzuola, J. Nodal line estimates for the second Dirichlet eigenfunction.  Journal of Spectral Theory (2019).  To appear.
  * [7] Behnke, H., and Goerisch, F. Inclusions for eigenvalues of selfadjoint problems.  In Topics in validated computations (Oldenburg, 1993), vol. 5 of Stud. Comput. Math. North-Holland, Amsterdam, 1994, pp. 277–322.
  * [8] Betcke, T. The generalized singular value decomposition and the method of particular solutions.  SIAM J. Sci. Comput. 30, 3 (2008), 1278–1295.
  * [9] Betcke, T., and Trefethen, L. N. Reviving the method of particular solutions.  SIAM Rev. 47, 3 (2005), 469–491.
  * [10] Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B. Julia: A fresh approach to numerical computing.  SIAM Review 59, 1 (2017), 65–98.
  * [11] Carstensen, C., and Gedicke, J. Guaranteed lower bounds for eigenvalues.  Math. Comp. 83, 290 (2014), 2605–2629.
  * [12] Courant, R., and Hilbert, D. Methods of mathematical physics. Vol. I.  Interscience Publishers, Inc., New York, N.Y., 1953.
  * [13] Dahne, J., and Salvy, B. Computation of tight enclosures for Laplacian eigenvalues.  SIAM J. Sci. Comput. 42, 5 (2020), A3210–A3232.
  * [14] Damascelli, L. On the nodal set of the second eigenfunction of the Laplacian in symmetric domains in $\mathbb{R}^{N}$.  Atti Accad. Naz. Lincei Cl. Sci. Fis. Mat. Natur. Rend. Lincei (9) Mat. Appl. 11, 3 (2000), 175–181 (2001).
  * [15] Fairweather, G., and Karageorghis, A. The method of fundamental solutions for elliptic boundary value problems.  vol. 9. 1998, pp. 69–95.  Numerical treatment of boundary integral equations.
  * [16] Fournais, S. r. The nodal surface of the second eigenfunction of the Laplacian in ${\bf R}^{D}$ can be closed.  J. Differential Equations 173, 1 (2001), 145–159.
  * [17] Fox, L., Henrici, P., and Moler, C. Approximations and bounds for eigenvalues of elliptic operators.  SIAM J. Numer. Anal. 4 (1967), 89–102.
  * [18] Freitas, P., and Krejčiřík, D. Unbounded planar domains whose second nodal line does not touch the boundary.  Math. Res. Lett. 14, 1 (2007), 107–111.
  * [19] Gershgorin, S. A. Über die Abgrenzung der Eigenwerte einer Matrix.  Bulletin de l’Académie des Sciences de l’URSS. Classe des sciences mathématiques et naturelles, 6 (1931), 749–754.
  * [20] Goerisch, F. Ein Stufenverfahren zur Berechnung von Eigenwertschranken.  In Numerical treatment of eigenvalue problems, Vol. 4 (Oberwolfach, 1986), vol. 83 of Internat. Schriftenreihe Numer. Math. Birkhäuser, Basel, 1987, pp. 104–114.
  * [21] Golberg, M. A., and Chen, C. S. The method of fundamental solutions for potential, Helmholtz and diffusion problems.  In Boundary integral methods: numerical and mathematical aspects, vol. 1 of Comput. Eng. WIT Press/Comput. Mech. Publ., Boston, MA, 1999, pp. 103–176.
  * [22] Gómez-Serrano, J. Computer-assisted proofs in PDE: a survey.  SeMA J. 76, 3 (2019), 459–484.
  * [23] Gómez-Serrano, J., and Orriols, G. Any three eigenvalues do not determine a triangle.  J. Differential Equations 275 (2021), 920–938.
  * [24] Gopal, A., and Trefethen, L. N. New Laplace and Helmholtz solvers.  Proc. Natl. Acad. Sci. USA 116, 21 (2019), 10223–10225.
  * [25] Gopal, A., and Trefethen, L. N. Solving Laplace Problems with Corner Singularities via Rational Functions.  SIAM J. Numer. Anal. 57, 5 (2019), 2074–2094.
  * [26] Grebenkov, D. S., and Nguyen, B.-T. Geometrical structure of Laplacian eigenfunctions.  SIAM Rev. 55, 4 (2013), 601–667.
  * [27] Grieser, D., and Jerison, D. Asymptotics of the first nodal line of a convex domain.  Invent. Math. 125, 2 (1996), 197–219.
  * [28] Hochman, A., Leviatan, Y., and White, J. K. On the use of rational-function fitting methods for the solution of 2D Laplace boundary-value problems.  J. Comput. Phys. 238 (2013), 337–358.
  * [29] Hoffmann-Ostenhof, M., Hoffmann-Ostenhof, T., and Nadirashvili, N. The nodal line of the second eigenfunction of the Laplacian in ${\bf R}^{2}$ can be closed.  Duke Math. J. 90, 3 (1997), 631–640.
  * [30] Jerison, D. The first nodal line of a convex planar domain.  Internat. Math. Res. Notices, 1 (1991), 1–5.
  * [31] Jerison, D. The diameter of the first nodal line of a convex domain.  Ann. of Math. (2) 141, 1 (1995), 1–33.
  * [32] Johansson, F. Arb: efficient arbitrary-precision midpoint-radius interval arithmetic.  IEEE Transactions on Computers 66 (2017), 1281–1292.
  * [33] Kato, T. On the upper and lower bounds of eigenvalues.  J. Phys. Soc. Japan 4 (1949), 334–339.
  * [34] Kennedy, J. B. Closed nodal surfaces for simply connected domains in higher dimensions.  Indiana Univ. Math. J. 62, 3 (2013), 785–798.
  * [35] Kennedy, J. B. A toy Neumann analogue of the nodal line conjecture.  Arch. Math. (Basel) 110, 3 (2018), 261–271.
  * [36] Lehmann, N. J. Optimale Eigenwerteinschliessungen.  Numer. Math. 5 (1963), 246–272.
  * [37] Lin, C. S. On the second eigenfunctions of the Laplacian in ${\bf R}^{2}$.  Comm. Math. Phys. 111, 2 (1987), 161–166.
  * [38] Liu, X. A framework of verified eigenvalue bounds for self-adjoint differential operators.  Appl. Math. Comput. 267 (2015), 341–355.
  * [39] Liu, X. Explicit eigenvalue bounds of differential operators defined by symmetric positive semi-definite bilinear forms.  J. Comput. Appl. Math. 371 (2020), 112666, 7.
  * [40] Liu, X., and Oishi, S. Verified eigenvalue evaluation for the Laplacian over polygonal domains of arbitrary shape.  SIAM J. Numer. Anal. 51, 3 (2013), 1634–1654.
  * [41] Melas, A. D. On the nodal line of the second eigenfunction of the Laplacian in ${\bf R}^{2}$.  J. Differential Geom. 35, 1 (1992), 255–263.
  * [42] Moler, C. B., and Payne, L. E. Bounds for eigenvalues and eigenvectors of symmetric operators.  SIAM J. Numer. Anal. 5 (1968), 64–70.
  * [43] Nakao, M. T., Plum, M., and Watanabe, Y. Numerical verification methods and computer-assisted proofs for partial differential equations, vol. 53 of Springer Series in Computational Mathematics.  Springer, Singapore, 2019.
  * [44] Payne, L. E. Isoperimetric inequalities and their applications.  SIAM Rev. 9 (1967), 453–488.
  * [45] Payne, L. E. On two conjectures in the fixed membrane eigenvalue problem.  Z. Angew. Math. Phys. 24 (1973), 721–729.
  * [46] Plum, M. Eigenvalue inclusions for second-order ordinary differential operators by a numerical homotopy method.  Z. Angew. Math. Phys. 41, 2 (1990), 205–226.
  * [47] Plum, M. Explicit $H_{2}$-estimates and pointwise bounds for solutions of second-order elliptic boundary value problems.  J. Math. Anal. Appl. 165, 1 (1992), 36–61.
  * [48] Pütter, R. On the nodal lines of second eigenfunctions of the fixed membrane problem.  Comment. Math. Helv. 65, 1 (1990), 96–103.
  * [49] Read, W. W., Sneddon, G. E., and Bode, L. A series method for the eigenvalues of the advection diffusion equation.  ANZIAM J. 45, (C) (2003/04), C773–C786.
  * [50] Shewchuk, J. R. Adaptive precision floating-point arithmetic and fast robust geometric predicates.  vol. 18. 1997, pp. 305–363.  ACM Symposium on Computational Geometry (Philadelphia, PA, 1996).
  * [51] Still, G. Computable bounds for eigenvalues and eigenfunctions of elliptic differential operators.  Numer. Math. 54, 2 (1988), 201–223.
  * [52] Tanaka, K. A posteriori verification for the sign-change structure of solutions of elliptic partial differential equations.  arXiv preprint arXiv:2001.03854 (2020).
  * [53] Tucker, W. Validated numerics.  Princeton University Press, Princeton, NJ, 2011.  A short introduction to rigorous computations.
  * [54] Weinstein, A., and Stenger, W. Methods of intermediate problems for eigenvalues.  Academic Press, New York-London, 1972.  Theory and ramifications, Mathematics in Science and Engineering, Vol. 89.
  * [55] Yang, D.-H., and Guo, B.-Z. On nodal line of the second eigenfunction of the Laplacian over concave domains in ${\mathbb{R}}^{2}$.  J. Syst. Sci. Complex. 26, 3 (2013), 483–488.
  * [56] Yau, S.-T. Open problems in geometry.  In Differential geometry: partial differential equations on manifolds (Los Angeles, CA, 1990), vol. 54 of Proc. Sympos. Pure Math. Amer. Math. Soc., Providence, RI, 1993, pp. 1–28.
  * [57] You, C., Xie, H., and Liu, X. Guaranteed eigenvalue bounds for the Steklov eigenvalue problem.  SIAM J. Numer. Anal. 57, 3 (2019), 1395–1410.

Joel Dahne
---
Department of Mathematics
Uppsala University
Lägerhyddsvägen 1, 752 37, Uppsala, Sweden
Email<EMAIL_ADDRESS>
Javier Gómez-Serrano
Department of Mathematics
Brown University
Kassar House, 151 Thayer St.
Providence, RI 02912, USA
and
Departament de Matemàtiques i Informàtica
Universitat de Barcelona
Gran Via de les Corts Catalanes, 585
08007, Barcelona, Spain
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
Kimberly Hou
Department of Mathematics
Princeton University
Fine Hall, Washington Rd,
Princeton, NJ 08544, USA
Email<EMAIL_ADDRESS>