The hot spots conjecture can be false: Some numerical examples

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.


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). ‡ Author to whom any correspondence should be addressed.
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 ∂ t u = ∆u for t → ∞ with homogeneous Neumann boundary condition ∂ ν 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.  The numerical solution of the heat equation ∂ t u = 1 10 ∆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 = const. Mathematically, we have to find a solution u = 0 and the smallest k ∈ R >0 such that the Helmholtz equation ∆u + k 2 u = 0 in D with ∂ ν u = 0 on the boundary Γ. All solutions k are called non-trivial interior Neumann eigenvalues and λ 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 = λ 0 < λ 1 ≤ λ 2 ≤ λ 3 ≤ . . . when D is a bounded planar domain with Lipschitz boundary (see for example [16, p. 449] and the references therein, specifically [33]). If λ 1 < λ 2 and u(0, x), u 1 = 0, then u(t, x) = e λ 1 t u(0, x), u 1 + 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 ⊂ R 2 be an open connected bounded domain with Lipschitz boundary Γ. Then: C1: For each eigenfunction u 2 (x) corresponding to λ 2 which is not identically zero, we have inf x∈Γ u 2 (x) < u 2 (y) < sup x∈Γ u 2 (x) ∀y ∈ D .
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 nonlinear 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.

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]).

Notations
We consider a bounded Lipschitz domain D ⊂ R 2 with one hole. The outer boundary Γ 1 is assumed to be sufficiently smooth that is oriented counter-clockwise and a sufficiently smooth inner boundary Γ 2 that is oriented clockwise. The normal ν 1 on the boundary Γ 1 is pointing into the unbounded exterior E. The normal ν 2 on the boundary Γ 2 is pointing into the bounded exterior I. We refer the reader to Figure 2. The boundary of D is given by Γ = Γ 1 ∪ Γ 2 (Γ 1 ∩ Γ 2 = ∅). Note that we also consider bounded domains without a hole. In this case, we have Γ 2 = ∅ and hence Γ = Γ 1 and I = ∅.

Boundary integral equation
The solution to the Helmholtz equation ∆u + k 2 u = 0 (reduced wave equation) in the domain D for a given wave number k with Im(k) ≥ 0 is given by (see [12,Theorem 2 which can be written in our notation as where Φ k (x, y) = iH (1) 0 (k x − y )/4, x = y denotes the fundamental solution of the Helmholtz equation in two dimensions (see [12, p. 66]). Here, H (1) 0 denotes the firstkind Hankel function of order zero. We denote u(y) for y ∈ Γ 1 as u 1 (y) and similarly u(y) for y ∈ Γ 2 as u 2 (y). Hence, we can write (1) as where we also used the homogeneous Neumann boundary conditions ∂ ν 1 u = 0 and ∂ ν 2 u = 0. We rewrite (2) as where we used the notation for the acoustic double layer potential with density ψ 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 ∈ D approach the boundary Γ 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 where we used the notation for the double layer operator (see [12, p. 41]). Here, Ω 1 (x) denotes the interior solid angle at a point x on Γ 1 . When the boundary is smooth at this point, then Ω 1 (x) = 1/2. In fact, Ω 1 (x) is 1/2 almost everywhere for Lipschitz domains.
Similarly, we obtain for x ∈ D approaching the boundary Γ 2 and using the jump relation for the double layer operator We can rewrite (4) and (5) as a 2 × 2 system of boundary integral equations in the form  where I and I B denotes the identity and the 2 × 2 block identity operator, respectively. Hence, we have to numerically solve the non-linear eigenvalue problem (6) written as 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 K(k) maps from denotes a Sobolev space of order s ∈ R on the domain Γ which are defined via Bessel potentials (see [28, pp. 75-76] for more details). The operator M(k) = C· I B + K(k) is Fredholm of index zero for k ∈ C\R ≤0 and therefore the theory of eigenvalue problems for holomorphic Fredholm operator-valued functions applies to M(k).

Discretization
In this section, we explain how to discretize (6) using quadratic interpolation of the boundary, but using piecewise quadratic interpolation with α = (1 − 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 where the matrix is of size 3· 2· n f ×3· 2· 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 First, we subdivide the boundary Γ 1 into n f pieces denoted by ∆ j with j = 1, . . . , n f . A subdivision into four pieces is shown in Figure 3 for the unit circle. Then equation (8) can be equivalently written as For each j there exists a unique map m j which maps from the standard interval σ = [0, 1] to ∆ j . Then, we can apply a simple change of variables to each integral over ∆ j giving with the Jacobian given by J(s) = ∂ 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 with u = 1 − s. Here, γ mod δ = γ − ⌊ γ δ ⌋· δ with ⌊· ⌋ the floor function. The nodes v ℓ (ℓ = 1, . . . , 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). An example how the approximation of ∆ 1 via a quadratic interpolation polynomial using two vertices and the midpoint looks like is shown in Figure 5. Next, we define the 'collocation nodes' v j,k by v j,k = m j (q k ) for j = 1, . . . , n f and for k = 1, 2, 3 where q 1 = α, q 2 = 1/2, and q 3 = 1 − α with 0 < α < 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 α the overall convergence rate can be improved. The first three collocation nodes on the approximated boundary for the unit circle using α = (1 − 3/5)/2 are shown in Figure 6.
The unknown function u 1 ( m j (s)) is now approximated on each of the j pieces by a quadratic interpolation polynomial of the form 3 k=1 u 1 ( m j (q k )) L k (s) which can be written as 3 k=1 u 1 ( v j,k ) L k (s) where the Lagrange basis functions are given by  with u = 1 − s. We obtain with r(x) the residue which is due to the different approximations. We set r( v i,ℓ ) = 0 and since Ω 1 ( v i,ℓ ) = 1/2 always by the choice of the collocation nodes, we obtain the linear system of size 3n f × 3n f with the resulting integrals which will be approximated by the adaptive Gauss-Kronrod quadrature (see [36]). This can be written abstractly as M(k) u = 0. Note that the integrand of the integral of the right-hand side (9) can easily by written down as the first-kind Hankel function of order one and with Note that the Jacobian cancels out. A word has to be spent on the following issue: When P = Q with P = v i,ℓ and Q = 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 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 smooth i,ℓ = 0. The integral I sing i,ℓ (a singularity is present) can be rewritten as a sum of integrals without any singularity. This is due to the fact that we have [35, p. 363]) and hence, we approximately use for all i, ℓ: (that is, the row sum of the matrix M(0) obtained from the discretization of the double layer for the Laplace equation shall be −1/2) and hence, we can compute I sing i,ℓ as 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 M(k). Finally, if we use a domain with a hole, we obtain the system of size m × m = 3· 2· n f ×3· 2· n f when including the boundary Γ 2 and the unknown function u 2 . Written abstractly we obtain the non-linear eigenvalue problem Of course, the extension to more than one hole is obvious. In this case, the matrix M(k) within (11) will be of size (q − 1)· 2· n f × (q − 1)· 2· n f where q denotes the number of holes within the domain.

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π-periodic contour γ 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 (µ, 0) (in order to exclude the eigenvalue zero, we choose µ > R). In this case, we have ϕ(t) = µ + R cos(t) + iR sin(t) which satisfies ϕ ∈ C ∞ . The number of eigenvalues including their multiplicity within the contour γ is denoted by n(γ). With the randomly chosen matrixV ∈ C m×ℓ with m ≫ ℓ ≥ n(γ) the two contour integrals of the form over the given contour γ are now rewritten as and approximated by the trapezoidal rule yielding where the parameter N is given and the equidistant nodes are t j = 2πj/N, j = 0, . . . , 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 A 0,N = VΣW H with V ∈ C m×ℓ , Σ ∈ C ℓ×ℓ , and W ∈ C ℓ×ℓ . Then, we perform a rank test for the matrix Σ = diag(σ 1 , σ 2 , . . . , σ ℓ ) for a given tolerance ǫ = tol rank (usually ǫ = 10 −4 ). That is, find n(γ) such that , and W 0 = (W ij ) 1≤i≤ℓ,1≤j≤n(γ)) and compute the n(γ) eigenvalues k i and eigenvectors s i of 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.

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 for an arbitrary point x ∈ 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].

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 threedimensional 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 < α < 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 α = (1 − 3/5)/2 (a Gauss-quadrature point within the interval [0, 1]), we have not investigated this any further.

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 (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 For the Beyn algorithm we use the parameters N = 24, R = 1/2, µ = 2, and ℓ = 10 for various number of faces n f and number of collocation points n c . With the definition of the absolute error E (i) n f of the i-th eigenvalue approximation, we compute the error of the first non-trivial eigenvalue E (1) n f of our method compared with (12). Additionally, we define the estimated error of convergence EOC (i) = log(E (2) of the i-th eigenvalue approximation and compute EOC (1) . As we can see in Table 1, the first nontrivial Neumann eigenvalue can be made accurate up to eleven digits. The estimated 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. 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 ′ 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 given by 3.831 705 970 207 512. Note that the first root of J ′ 0 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, µ = 3, and ℓ = 10.
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.
. In general, we use a resolution of 100 × 100 equidistantly distributed points, here within 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 nontrivial interior Neumann eigenvalue as a red and blue dot, respectively. 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. The boundary of the ellipse is given in parametric form as (6 cos(t)/5, sin(t)) ⊤ with t ∈ [0, 2π). We use the parameters as before for Beyn's algorithm except µ = 1.5 and consider 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) + ε cos(2t), sin(t)) ⊤ with t ∈ [0, 2π), where the parameter ε 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 ε = 0.1 and ε = 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 π and 4π/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.
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.   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 π and 4π/3.

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 ε = 0.3, the peanut-shaped domain and the apple-shaped domain. The boundaries of the last two domains are given parametrically as cos 2 (t) + sin 2 (t)/4 (cos(t), sin(t)) ⊤ and 0.5+0.4 cos(t)+0.1 sin(2t) 1+0.7 cos(t) (cos(t), sin(t)) ⊤ with t ∈ [0, 2π), respectively (see [41] for its use). We use the parameters as before with µ = 1.5 for the first two domains and µ = 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 1.3 .
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  All three eigenfunctions including the location of the maximal and minimal value are shown in Figure 11. 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 = [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 √ 1.475 621 845 ≈ 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. In Figure 12 we show the corresponding eigenfunction. We also tried two different domains 2 and 3 as shown in Figure 4 b) and c). Using the same parameters as for the L-shaped domain , we obtain the first non-trivial Neumann eigenvalues 2.725 559 and 2.207 276, respectively. Interestingly, the    1.0716, 1.1645, 1.1745, 1.2202, 1.3329, and 1.6816 for 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.

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  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 µ = 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. In Figure 13 we show the first three eigenfunctions for 2.1 corresponding to the nontrivial interior Neumann eigenvalues 0.822 253, 1.504 648, and 2.096 773. To compute the third eigenvalue which has multiplicity two, we used µ = 2. For the first eigenfunction plot, we also added the extreme values.
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   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: Let E(x, y, a, b, t 1 , t 2 ) be the ellipse centered at (x, y) ⊤ with semi-axis a and  Figure 15. A more advance bounded domain with a hole inspired by the work of Burdzy [9, Figure 1].
b constructed for φ ∈ [t 1 , t 2 ) using the parametrization (x + a cos(φ), y + b sin(φ)) ⊤ . Note that we also allow t 1 > t 2 to guarantee the needed orientation of the curve.
We use the parameter µ = 0.8, ℓ = 20 and 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. 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 1.6 . As we can see again in Figure 17, the values inside the bump  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. As we can see in Figure 18, the location of the hot spots of the eigenfunction within 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 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, −π) and E(0, 3/2, 1, 1/2, 0, π) and its mirror counterpart smaller. We use E(0, 4, 1, 1, 0, −π) and E(0, 3/2, 1, 1, 0, π) instead to obtain C 4 . We obtain the first non-trivial interior Neumann eigenvalue 0.384 715 and its corresponding eigenfunction within 1.8 shown in Figure  19. 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 δ > 0. The first half of the outer boundary is given by the pieces E(4, 3 + δ, 1, δ, −π/2, −π), E(2, 3 + δ, 1, δ, 0, π), E(0, 3 + δ, 1, δ, 0, −π), E(−2, 3 + δ, 1, δ, 0, π), E(−4, 3 + δ, 1, δ, 0, −π/2), and E(−4, 0, 3, 3, π/2, 3π/2). Rotating this half by π yields the second half of the outer boundary. The first half of the inner boundary is given by the pieces E(4, 5/2 − δ, 1, δ, π/2, π), E(2, 5/2 − δ, 1, δ, 0, −π), E(0, 5/2 − δ, 1, δ, 0, π), E(−2, 5/2 − δ, 1, δ, 0, −π), E(−4, 5/2 − δ, 1, δ, 0, π/2), and E(−4, 0, 5/2, 5/2, π/2, 3π/2). Rotating this half by π 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 δ = 1 yields the domain C 4 .
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 1.8 . The corresponding non-trivial interior Neumann eigenvalue is 0.563 329. 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 1.8 . 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 alsoC 4 andĈ 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)· (sin(t), cos(t)) ⊤ with r i (t) = a i (1 + 1/2 sin(ωt· I [(ω/2−1)π/ω,(ω/2+1)π/ω]∪[(3ω/2−1)π/ω,(3ω/2+1)π/ω] )), t ∈ [0, 2π) with a 1 = 1 for the outer boundary and a 2 = 0.9 for the inner boundary. Here I denotes the indicator function. We use ω = 4 and ω = 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 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. 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 ℵ max and ℵ 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 ℵ max > 1 and/or ℵ 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) ⊤ and (0, −0.6) ⊤ , respectively. For the next three domains, we use (0, 0.69) ⊤ and (0, −0.69) ⊤ and for the last domain, we use (0, 0.6875) ⊤ and (0, −0.6875) ⊤ .

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) ⊤ with radius R, yields a domain with two holes, say C 1,R . Using R = 0.1, R = 0.15, and 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 . 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,R and mirror the circular hole at the y-axis. This yields the domain C 1,R . Using the same parameters as before yields the results shown in Figure 26. As we can see now, the maximum and minimum remain inside all the three domains C 1,0.1 , C 1,0.15 , and 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 mir 1,0.1 , C mir 1,0.15 and C mir 1,0.2 , respectively. The results are presented in Figure 27. As expected, we obtain the minimal and maximal value inside the domain for the two domains C mir 1,0.1 and C mir 1,0.15 whereas the maximum is on the boundary for the domain C mir 1,0.2 . The domains with five holes C mir 1,0.1 , C mir 1,0.15 and C mir 1,0.2 are constructed by mirroring the domains C 1,0.1 , C 1,0.15 and C 1,0.2 at the x-axis. The results are shown in Figure 28. 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.

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.