Computation of Laplacian eigenvalues of two-dimensional shapes with dihedral symmetry

We numerically compute the lowest Laplacian eigenvalues of several two-dimensional shapes with dihedral symmetry at arbitrary precision arithmetic. Our approach is based on the method of particular solutions with domain decomposition. We are particularly interested in asymptotic expansions of the eigenvalues $\lambda(n)$ of shapes with $n$ edges that are of the form $\lambda(n) \sim x\sum_{k=0}^{\infty} \frac{C_k(x)}{n^k}$ where $x$ is the limiting eigenvalue for $n\rightarrow \infty$. Expansions of this form have previously only been known for regular polygons with Dirichlet boundary condition and (quite surprisingly) involve Riemann zeta values and single-valued multiple zeta values, which makes them interesting to study. We provide numerical evidence for closed-form expressions of higher order $C_k(x)$ and give more examples of shapes for which such expansions are possible (including regular polygons with Neumann boundary condition, regular star polygons and star shapes with sinusoidal boundary).


Introduction
Let Ψ be a function defined on a domain Ω ⊂ R 2 that satisfies the Laplace eigenvalue equation where ∆ = ∂ 2 x +∂ 2 y denotes the Laplacian in two-dimensional Euclidean space and λ ∈ R. We consider the case when the two-dimensional shape Ω has the symmetry group of a regular n-gon, and Ψ has either Dirichlet or Neumann boundary condition on ∂Ω Ψ| ∂Ω = 0 (Dirichlet) , ∂ ⃗ n Ψ| ∂Ω = 0 (Neumann) , where ∂ ⃗ n denotes the normal derivative. By combining Eq. (1.1) with the boundary condition, one obtains a discrete spectrum of eigenvalues 0 < λ 1 ≤ λ 2 ≤ λ 3 ≤ . . . , (in the case of Neumann boundary condition we additionally have λ 0 = 0). The eigenvalues are invariant under translations and rotations but do depend on the area of the considered shape. More precisely, if Ω ′ is obtained from Ω by a homothety, then one has where A denotes the area of the corresponding shape. Keeping the area constant is therefore crucial for investigating 1/n expansions. Throughout this paper, we will always normalize the domains to have area π, which means that the shapes will approach the unit disk in the limit as n → ∞. We start by considering regular polygons with Dirichlet boundary condition. Additional examples will be introduced in Section 7. There are so far only two regular polygons whose eigenvalues are known explicitly. These are the regular triangle where λ 1 = 4π/ √ 3 and the square with λ 1 = 2π, see, e.g., Pólya and Szegö [32] (one should however mention that some eigenmodes of the regular hexagon can be obtained by piecing together eigenmodes of regular triangles). The remaining regular polygons offer challenges, both analytically and numerically, due to the presence of non-analytic vertices (i.e. vertices that are not of the form π/n). It has been known from works of several authors [2,7,18,25,30] that the lowest eigenvalue of a regular n-gon with Dirichlet boundary condition can be approximated by a series in 1/n where x = j 2 0,1 is the liming eigenvalue (i.e., the lowest eigenvalue of the unit disk, which is given by the square of the first root of the Bessel function J 0 (x)). Since regular polygons are approaching the circle in the limit, it is natural that C 0 = 1. Interestingly, the remaining coefficients seemed to be expressible in closed-form as polynomials in x = j 2 0,1 of degree ⌊ i−3 2 ⌋ with coefficients expressible in terms of Riemann zeta values C 0 (x) = 1 , The first coefficients up to third order have been discovered by Grinfeld and Strang [18] in 2012 through deforming the circle to a regular polygon and then applying a technique called calculus of moving surfaces. Three years later, Boady in his PhD thesis [7] managed to compute two more terms also using methods from the calculus of moving surfaces. The seventh and eighth coefficients were recently found by the second author [25] using numerical methods similar to the ones presented in this report.
We managed to obtain eight higher order coefficients in this project C 9 (x) = ζ(9)( 9 4 x 3 + 104x 2 + 438x − 1020) + ζ(3) 3 (240x + 96) , C 10 (x) = ζ(7)ζ(3)(x 3 + 39x 2 − 24x + 144) + ζ(5) 2 (x 3 − 6x 2 − 12x + 72) , . (The expressions of the coefficients C 14 (x), C 15 (x), and C 16 (x) are given in the Appendix, the numerical expressions as well as the underlying eigenvalue data are provided as an attachment to this paper). In order to determine these expressions we computed the eigenvalues of 650 n-gons to at least 980 digits precision. Through this we also discovered the appearance of single-valued multiple zeta values (MZVs) in the expansion coefficients, Figure 1. Fundamental region of regular polygons for fully symmetric eigenfunctions starting at eleventh order. In case of C 16 (x) the basis of single-valued multiple zeta values was found with help of the program HyperlogProcedures developed by Oliver Schnetz. We provide a brief introduction to MZVs and the definitions of the single-valued MZVs that appear in our expressions in the Appendix. We should note that the approach used in this paper does not produce explicitly proven results but rather gives conjectures with very high numerical evidence. We did however manage to prove some of the results for regular polygons with Dirichlet boundary condition in the companion paper [2]. Namely, we explicitly derived the results up to (and including) C 14 and proved that that the expansion coefficients are expressions over the space of multiple zeta values. However, the expressions for C 15 and C 16 can only be considered conjectural at this point.

The Method of Particular Solutions
Because the eigenfunctions considered in this work are dihedrally symmetric, it is sufficient to compute them inside a triangular subdomain (which we refer to as fundamental domain) instead of working with the full shape (see Fig. 1 for an illustration of the fundamental domain for a regular polygon). A numerical method for computing eigenvalues of triangular shapes is given by the method of particular solutions (MPS). The MPS has been introduced by Conway [10] in 1961 and established by a famous paper of Fox, Henrici and Moler [16] in 1966 on computing eigenvalues of L-shape domains.
The main idea of the MPS is to expand Ψ for the spectral parameter k = √ λ as a Fourier-Bessel series (in polar coordinates) and J m denotes the Bessel function (see for example [35]). The basis functions ψ ν have the useful property that they can serve as eigenfunctions along an unbounded wedge by a wellmade choice of m ν . Consider for example a wedge with angle α with the vertex at the origin. Dirichlet boundary conditions impose that the function has to vanish along the boundary. If one chooses m ν = νπ/α for the odd eigenfunction basis, one obtains functions that trivially fulfill the boundary condition along both edges of the wedge. This property (combined with the exponential decay of the basis functions) makes the MPS very useful when computing eigenvalues of triangular shapes: the boundary condition on two of the three edges can be trivially fulfilled by the construction of ψ ν . To obtain the spectral parameter k one truncates the expansion of Ψ to some finite order N and searches for the parameter k for which the expansion vanishes on a discrete set of points on the remaining edge up to the desired numerical precision (this approach is often referred to as point-matching). This results in a linear system of equations which is square if the number of point-matching points is also chosen to be N . The parameter k now corresponds to an approximation of the true spectral parameter if the coefficients are essentially invariant under the choice of matching points. In view of this, one computes the value of k for which the determinant becomes zero (i.e., when the matrix becomes singular) using root-finding algorithms (we used the secant method) [16]. By increasing N , one can obtain better approximations and hence better precision for the eigenvalue. This approach has been applied by the second author [24,25] who expanded from the vertex with angle β/2 (see Fig. 1) and applied point-matching to the edge ∂Ω 2 to compute eigenvalues of polygons with n ≤ 150 to 50 digits of precision. We remark that the MPS has the limitation that the point-matching matrix becomes ill-conditioned for large N . This has been further analyzed by Betcke and Trefethen [4,5] in 2005 who proposed an improved approach that is referred to as modified MPS. Their main idea has been to compute an additional matrix that consists of (randomly chosen) points in the interior of the considered region. This matrix is designed to ensure that the eigenfunction is non-zero in the interior. By minimizing the lowest generalized singular value of these two matrices, one can reliably locate the eigenvalues of a given shape because spurious solutions are excluded and N can be chosen as required, without conditioning the matrix too poorly. The modified MPS might be regarded as a successor of the MPS and is superior in most applications. For the computation of eigenvalues of relatively simple shapes at arbitrary precision arithmetic, the original MPS can however still prove itself to be very competitive because the ill-conditioning can be overcome by increasing the working precision and because the generalized singular value decomposition is very computationally expensive compared to the LU-decomposition that is required for the determinant computation. Additionally, the function of the determinant w.r.t. k becomes almost linear close to the eigenvalues which makes the location of the roots very efficient [24].
We further remark that recent progress has been achieved in making the results of the MPS rigorous (i.e., with certified error bounds) [11,20]. For our application, certifying the eigenvalues yields no benefit because the derivation of the closed-form expressions using the LLL algorithm is non-rigorous (unless one knows a bound of the height of the coefficients in advance, which we are unaware of). Additionally, specifically certifying the eigenvalues results in a loss of precision.

Domain Decomposition
Despite the absolute convergence of the MPS (i.e., one always gets better estimates by increasing N as long as the working precision is sufficiently increased) the algorithm which has been applied by the second author in [24,25] has the disadvantage that the convergence rate (which is the eigenvalue precision per amount of matching points) drastically decreases for polygons with more edges (and hence "thinner" fundamental regions). To overcome  [13] in 1983 who decomposed domains into several subdomains Ω k with eigenfunctions f k and minimized the quotient of the functionals where Γ kl denote straight intersecting boundaries. This method has been further developed by Driscoll [14] in 1997 to compute eigenmodes of the famous isospectral drums of Gordon, Webb and Wolpert [17]. Betcke [3] in 2006 formulated the method of Descloux and Tolley as a GSVD problem which does not require the evaluation of the boundary-and domain integrals.
Our approach is based on the idea to treat the intersection of regions as an additional point-matching condition. We chose to decompose the fundamental region of regular polygons into two regions. Then, instead of minimizing Eq. (3.1), we explicitly matched the eigenfunctions of both regions and their derivatives on a discrete set of points along the intersecting line. We do not claim this approach to be new (in fact it was already mentioned in Driscoll's paper [14]) however we are unaware of any author proving the success of this method.

Outline of the Algorithm
As shown in Fig. 2, the fundamental region can be decomposed into two regions I & II. For the first region one chooses the eigenfunctions to be and expands of the vertex with angle α/2. To fulfill the Neumann boundary condition at angle θ = 0 (which is required for fully symmetric eigenmodes), the angular derivative of the wave function (which is a sine with a factor that does not matter here) has to vanish. This is always fulfilled since sin(m ν · 0) = 0 for all m ν . The second Neumann boundary condition is fulfilled if sin(m ν α 2 ) = 0 from which follows Similarly for the second domain one expands from the β/2-vertex (and places it at the origin). We denote this coordinate system with a tilde. The space of eigenfunctions now becomes The point-matching matrix now has to treat three different boundary conditions: (1) The wave functions of both regions have to match along the intersecting line (2) The normal derivatives of both wave functions have to match along the intersecting line . . .
where the minus sign appears because the expansions have opposite orientation along the intersection line. We computed the derivatives explicitly by differentiating the Bessel and sine functions.   All these conditions can then be combined into a point-matching matrix We chose the height of the intersecting line to be equal to the length of ∂Ω 3 so that region II becomes close to a square. We have also experimented with different ways of decomposition (such as choosing the intersection line crosswise through the region, which eliminates the third point-matching condition), but none of them provided better results than the presented approach. In order to overcome the ill-conditioning of the point-matching matrix we performed the computations with about 1.2N digits working precision where N is the size of the point-matching matrix (this choice was obtained through trial and error, one could also attempt to choose N based on the asymptotic decay of the series, which does however require a growth condition for c ν ). The points along all boundaries (exterior and interior) were distributed using the Chebyshev distribution. By empirical testing we chose the number of matching points for all three boundary conditions equally to be N/3 each while the first eigenfunction is expanded up to N I = N/4 and N II = 3N/4. The advantage of this approach becomes clear in Fig. 3: the algorithm with domain decomposition is almost unaffected by the number of polygon edges n while the original point-matching algorithm without domain decomposition suffers from a heavy decrease of the convergence rate for increasing n. This allowed us to compute eigenvalues for hundreds of regular polygons to high precision.

Details on the numerical implementation
We implemented the algorithm in the Julia programming language [6]. We made use of Arb [21] and its Julia interface implemented in Nemo [15]. Arb is a highly optimized C library for arbitrary precision arithmetic. We used Arb to evaluate the occurring special functions. Especially the evaluation of Bessel functions is computationally demanding (in fact, most computational time was spent on computing these functions) which makes an efficient implementation very important (see [22] for some benchmarks of Arb's Bessel implementation). To compute the determinants we used a wrapper to the function arb mat approx lu of Arb which also runs considerably faster than the default Julia determinant algorithm [23]. In general, Arb is designed for interval-arithmetic (which means that every float has an error-bound attached so that rounding errors during arithmetic operations are taken into account). However, for our implementation interval-arithmetic is counterproductive because the point-matching matrix is ill-conditioned and empirical tests have revealed that interval-arithmetic might cancel significant digits too pessimistically during the determinant computation (as it ensures that the result is correct up to its displayed precision). It is therefore important to note that arb mat approx lu is one of the few functions of Arb where approximate (non-interval) arithmetic is used. The computing cluster at which the computation were performed consisted of nodes with 2x Intel Xeon E5-2680 v4 @ 2.40GHz CPUs with 14 (physical) cores each and 128GB of RAM per node. In total, 39 nodes were available. The computations were queued using HTCondor [29] which offers priority-based distribution of hardware resources between multiple users. This means that the hardware was shared with other users which might have increased CPU-times. Each time the determinant of a point-matching matrix has been computed (within the secantmethod) the result (and other significant parameters) were stored and the job re-queued. This approach allowed other users to occupy resources as well and the machines to perform required reboots. If however no other user was using the computing cluster, all resources could be used to prevent idle time. The largest computed matrices were of size 3039 × 3039 with 3647 digits precision per entry. The computation of these matrices (and their determinant) required around 45 hours of wall-clock time running on 9 threads (in our latest version we computed the point-matching matrix in parallel and then computed the determinant using Arb). To store such a matrix and compute its determinant, about 20 GB of RAM was needed. The computations took several months and required 2.7 million threadhours in total (the amount of physical core-hours without hyper-threading is impossible to determine).

Derivation of the coefficients
To compute the coefficients of the 1/n-expansion, we used Richardson-extrapolation as it is described in [1]. The computed eigenvalues λ(n) for n-sided polygons can be expanded as a series λ(n) = C 0 + C 1 · n −1 + C 2 · n −2 + · · · + C N · n −N λ(n + 1) = C 0 + C 1 · (n + 1) −1 + C 2 · (n + 1) −2 + · · · + C N · (n + 1) −N . . . λ(n + N ) = C 0 + C 1 · (n + N ) −1 + C 2 · (n + N ) −2 + · · · + C N · (n + N ) −N were n is the lowest computed polygon and N is the total amount of computed polygons. Known coefficients were subtracted to increase the precision. This creates a square system of equations which we solved with LU-decomposition to obtain numerical estimates for the coefficients C i . Using this approach gives higher precision than a standard polynomial fit as it was applied in [25]. Still, we were unable to extract the coefficients with nearly the same precision as the computed eigenvalues. For example for the seventeenth coefficient (for which we could not find a closed-form solution) we get a precision of approximately 580 digits while the eigenvalues were computed to at least 980 digits. We tried multiple ways to extract the coefficients with higher precision (for example by trying to exploit the fact that the Richardson matrix is a Vandermonde-matrix (see [12]) or by performing the computations with rational arithmetic to not implement any rounding errors) but none of them yielded in better results.
An empirical observation is that one seems to get the best results by using about 0.65D eigenvalues for the Richardson extrapolation, where D denotes the (estimated) precision of the eigenvalues in digits. This is the reason why we have computed 650 polygons to about 1000 digits precision.
After computing the numerical expressions of the coefficients, we used Pari [34] to find closed-form guesses using the LLL-algorithm [28]. For regular polygons with Dirichlet boundary condition, we have also made use of the results of [2], which give closed forms for the coefficients of x 0 and x 1 , to cut down on the search space.

Regular Polygons with Neumann Boundary Condition.
The presented algorithm works equally well for regular polygons with Neumann boundary condition (one only has to slightly adapt the function space for the second region). We have computed the lowest Laplacian eigenvalue of 455 polygons with Neumann boundary condition to at least 700 digits precision. This computed data provides strong evidence that there is a 1/n-expansion for the Neumann case with coefficients (here x = j 2 1,1 ) C 0 (x) = 1 , Our results thus indicate that the appearance of single-valued MZVs is not restricted to the Dirichlet eigenvalues (although the results for Neumann boundary condition have not been formally proven yet).

Regular Star Polygons.
In order to study whether the expandability of eigenvalues also applies for other polygonal shapes with dihedral symmetry, we have also considered regular star polygons, which are non-convex regular polygons. For convenient labeling, we make use of the Schläfli notation: Any star polygon can be labeled as {n, q} where n corresponds to the number of vertices and q ≥ 2 is referred to as the density (see Fig. 4). The numbers q and n are also required to be relatively prime. The fundamental regions are then given by the triangles with angles ( π n , π(n−2q) 2n , π(n+2q−2) 2n ). One example of a regular star polygon is the pentagram which is labeled as {5, 2}. The lowest 79 eigenvalues of this shape have been computed by the second author [24] to 20 digits of precision. Despite that, we are unaware of any efforts to compute the eigenvalues of these shapes to significantly high precision. The reason for that is probably that they are relatively difficult to compute at high precision since they contain two non-analytic vertices. In fact, we had to compute expansions using three vertices as shown in Fig. 5. We have computed expansions of the eigenfunctions in each region to roughly equal orders: N I ≈ N II ≈ N III ≈ N/3 (we also chose equal numbers of points along each intersection line). The working precision for this triple-expansion approach was chosen to be 0.7N + 50.
We have computed the eigenvalues of star shapes in the case q = 2 with 12 ≤ n ≤ 263 to almost 300 digits precision and in the cases 3 ≤ q ≤ 8 to about 100 digits of precision. We find that the expansion coefficients of regular star shapes also involve zeta values and can be written as polynomials in q. For the first nine coefficients we conjecture that . Here the number A 7 is defined as . Note that these expressions also reproduce the expansion coefficients of regular convex polygons which correspond to q = 1. What seems particularly interesting is that the coefficients C 7 and C 8 involve alternating multiple zeta values which did not appear in the case of regular polygons.

Smooth Star Shapes.
We also investigated the asymptotic expansion of the eigenvalues of shapes with curved vertices which we call smooth star shapes. Smooth star shapes are cycloids with radius where n corresponds to the number of arcs, d is the height of an arc and R is the radius of the underlying circle (see Fig. 6). The eigenvalues of these shapes have previously been investigated by Laura in 1964 [27] and Wagner in 1971 [36] who derived the approximate formula with c = j 0,1 . This formula gives an approximation of the eigenvalue at least up to the first digit for d = 0.3. Smooth star shapes are relatively easy to construct and have been used to test numerical algorithms for the computation of Laplacian eigenvalues, see for example [9,19]. In our case, we investigate the behavior of smooth star shape domains with d = (1/n) m where m ∈ N. This choice of radii ensures that the considered shapes have the eigenvalue of the circle as their limiting eigenvalue. To compute the eigenvalues we simply used an expansion from the interior vertex. This might not be the most efficient approach (since the point-matching matrix becomes extremely ill-conditioned) however it was sufficient to investigate the first terms of the asymptotic expansion. We find that the eigenvalues of these (artificially constructed) shapes can also be written as a 1/n series where the expansion coefficients are now rational polynomials of x = j 2 0,1 . For example for m = 2 one obtains the rational polynomials C 0 (x) = 1 , (The remaining coefficients up to C 23 as well as the coefficients for m = 3 and m = 4 can be found in Appendix.) We omit the case m = 1 simply for convenience (this case is slightly more difficult to compute to high precision) as we do not expect any substantially different behavior.

Hypocycloids.
A hypocycloid is the locus of a point on a small circle rotating along the interior of the boundary of a larger circle (see Fig. 7). In case when the ratio between the radii of the larger and smaller circles is an integer n, the resulting hypocycloid has the n-sided dihedral symmetry. Interestingly, these shapes are related to group theory: the traces of all matrices in SU(n) lie inside an n-sided hypocycloid as was shown by Kaiser [26]. Hypocycloids can be parameterized by where r = R/n is the radius of the smaller (rotating) circle. In the limit as n → ∞ one thus obtains the parametric equations of a circle of radius R. It seems therefore natural to investigate if the eigenvalues of hypocycloids can be expressed as a 1/n-expansion, similar to what was shown in the previous sections. The vertices of hypocycloids are called cusps, i.e., the internal angle at each vertex is zero. To compute the eigenvalues of a hypocycloid one can compute the expansions from the interior angle and from the cusp and match the two eigenfunctions at the overlap (see Fig. 8). Using the Fourier-Bessel expansion ansatz it is not possible to directly fulfill the boundary condition along the curved arcs. The choice of the second eigenfunction (centered at the origin) is therefore the general Fourier-Bessel The conditions are now that the two eigenfunctions (including their derivatives) have to match along the intersection line and additionally that the second function has to be zero along the arc. As far as we are aware, so far no other author has applied the method of particular solutions using expansions from a cusp vertex. The MPS (combined with domain decomposition) however works decently and achieves convergence with around 1 digit of eigenvalue precision per 25 additional matching points. We have computed each hypocycloid eigenvalues from 5 to 100 edges for one week on a single thread, which resulted in around 65 digits precision per eigenvalue. If one fits the eigenvalues on a 1/n-polynomial-expansion, the expansion coefficients are converging towards constant values which are given by We have so far been unable to find closed-form expressions for these numbers. Using the symmetry of the expansion formulae for higher fully symmetric eigenvalues (see for example [31]), we found indications that the expansion coefficients C i are polynomials of degree in the limiting eigenvalue x m = j 2 0,m . Namely, by computing the first four fully symmetric eigenvalues of several hypocycloids, we found evidence that However, it remains unclear whether these numbers can be expressed in closed form. is known to converge for all s ∈ C that satisfy Re(s) > 1. The Riemann zeta function is arguably the most famous function in all of mathematics due to its tight connection to the distribution of prime numbers. For our treatment we restrict ourselves to the special values of ζ(s) at natural numbers s ∈ N, s > 1, which are sometimes simply called zeta values. For even zeta values, Euler has famously showed that where B n are Bernoulli numbers. For odd zeta values much less is known. It has been proven by Apéry in 1979 that ζ(3) is irrational. His proof has not yet been extended to any other odd zeta values. It is known that at least one of ζ(5), ζ(7), ζ(9), ζ(11) has to be irrational and that there exists an infinite number of irrational odd zeta values. However, it is still unknown whether all odd zeta values are irrational and whether π, ζ(3), ζ(5), . . . , ζ(2n + 1), . . . have any algebraic relations among them (altough it is widely believed that there are none).
As we have seen in Section 1, the arguments of products of zeta values that appear in the expansion coefficients of polygon eigenvalues add up to the index of the corresponding asymptotic coefficient. For example, the tenth coefficient for the Dirichlet case contains the products ζ(7)ζ(3) and ζ(5) 2 , and we have 7 + 3 = 2 · 5 = 10. This number is called the weight of a zeta product. Note that the product of two zeta values can be written as We can split the last sum into three terms One of the terms is given as a zeta value The remaining "crossing" terms are the so-called multiple zeta values ζ(a, b) and ζ(b, a). More generally, the multiple zeta values (MZVs) are defined as The decomposition is sometimes called the Nielsen reflection formula. Multiple zeta values satisfy a wide variety of interesting relations. For instance, the Q-linear span of MZVs inside R is closed under products, i.e., forms an algebra. An important subalgebra of MZVs is given by single-valued MZVs that were defined by Brown in [8]. These have been found to appear in amplitudes of Feynman diagrams of string theory. The relation between single-valued MZVs and odd zeta values is given by (9.8) ζ sv (2n + 1) = 2ζ(2n + 1) . + ( 3100857 2800 ζ(13)ζ(3) + 139433 80 ζ(11)ζ(5) + 11185 72 ζ(9)ζ(7) − 20ζ (7)