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

David Berghaus, Robert Stephen Jones, Hartmut Monien, Danylo Radchenko

► To cite this version:

David Berghaus, Robert Stephen Jones, Hartmut Monien, Danylo Radchenko. Computation of Laplacian eigenvalues of two-dimensional shapes with dihedral symmetry. Advances in Computational Mathematics, 2024, 50 (3), pp.38. 10.1007/s10444-024-10138-3. hal-04789060

HAL Id: hal-04789060

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

Submitted on 18 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Computation of Laplacian eigenvalues of two-dimensional shapes with dihedral symmetry

David Berghaus<fnref n="1" /> ID · Robert Stephen Jones<fnref n="2" /> · Hartmut Monien<fnref n="1" /> · Danylo Radchenko<fnref n="3" />

Received: 6 December 2022 / Accepted: 4 April 2024 / Published online: 2 May 2024
© The Author(s) 2024

## Abstract

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 \to \infty$. Expansions of this form have previously only been known for regular polygons with Dirichlet boundary conditions 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).

Communicated by: Tobin Driscoll

<figure><img src="image_2.png" /></figure>

David Berghaus
davidberghaus7@gmail.com

Robert Stephen Jones
rsjones7@yahoo.com

Hartmut Monien
hmonien@uni-bonn.de

Danylo Radchenko
danradchenko@gmail.com

<footnote n="1">1. Bethe Center, University of Bonn, Nussallee 12, Bonn 53115, Germany</footnote>

<footnote n="2">2. Independent Researcher, Sunbury 43074, OH, USA</footnote>

<footnote n="3">3. Laboratoire Paul Painlevé, University of Lille, Villeneuve d'Ascq F-59655, France</footnote>
**Keywords** Laplacian eigenvalues · Method of particular solutions · Regular polygons · Star shapes · Multiple zeta values

**Mathematics Subject Classification (2010)** 11M32 · 65N25 · 65N35

# 1 Introduction

Let $\Psi$ be a function defined on a domain $\Omega \subset \mathbb{R}^2$ that satisfies the Laplace eigenvalue equation

$$ - \Delta \Psi(x, y) = \lambda \cdot \Psi(x, y), \qquad (1) $$

where $\Delta = \partial_x^2 + \partial_y^2$ denotes the Laplacian in two-dimensional Euclidean space and $\lambda \in \mathbb{R}$. We consider the case when the two-dimensional shape $\Omega$ has the symmetry group of a regular $n$-gon, and $\Psi$ has either Dirichlet or Neumann boundary condition on $\partial\Omega$

$$ \begin{aligned} \Psi|_{\partial\Omega} &= 0 \quad (\text{Dirichlet}), \\ \partial_{\vec{n}}\Psi|_{\partial\Omega} &= 0 \quad (\text{Neumann}), \end{aligned} $$

where $\partial_{\vec{n}}$ denotes the normal derivative. By combining Eq. (1) with the boundary condition, one obtains a discrete spectrum of eigenvalues

$$ 0 < \lambda_1 \le \lambda_2 \le \lambda_3 \le \dots, $$

(in the case of Neumann boundary condition, we additionally have $\lambda_0 = 0$). The eigenvalues are invariant under translations and rotations but do depend on the area of the considered shape. More precisely, if $\Omega'$ is obtained from $\Omega$ by a homothety, then one has

$$ \lambda_i(\Omega)A(\Omega) = \lambda_i(\Omega')A(\Omega'), \qquad (2) $$

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 $\pi$, which means that the shapes will approach the unit disk in the limit as $n \to \infty$.

We start by considering regular polygons with Dirichlet boundary conditions. Additional examples will be introduced in Sect. 8. There are so far only two regular polygons whose eigenvalues are known explicitly. These are the regular triangle where $\lambda_1 = 4\pi/\sqrt{3}$ and the square with $\lambda_1 = 2\pi$, see, e.g., Pólya and Szegő [1] (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 $\pi/n$). It has been known from works of several authors [2–6] that the lowest eigenvalue of a regular $n$-gon with
Dirichlet boundary condition can be approximated by a series in $1/n$

$$ \lambda_1(n) \sim x \sum_{k=0}^{\infty} \frac{C_k(x)}{n^k}, \qquad (3) $$

where $x = j_{0,1}^2$ is the limiting 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_{0,1}^2$ of degree $\lfloor \frac{i-3}{2} \rfloor$ with coefficients expressible in terms of Riemann zeta values

$$ \begin{aligned} C_0(x) &= 1, \\ C_1(x) &= 0, \\ C_2(x) &= 0, \\ C_3(x) &= 4\zeta(3), \\ C_4(x) &= 0, \\ C_5(x) &= \zeta(5)(-2x + 12), \\ C_6(x) &= \zeta(3)^2(4x + 8), \\ C_7(x) &= \zeta(7)(-\frac{1}{2}x^2 - 12x + 36), \\ C_8(x) &= \zeta(3)\zeta(5)(2x^2 + 8x + 48). \end{aligned} $$

The first coefficients up to third order have been discovered by Grinfeld and Strang [3] 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 [4] 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 [5] using numerical methods similar to the ones presented in this report.

We managed to obtain eight (conjectural) closed-form expressions of higher-order coefficients in this project

$$ \begin{aligned} C_9(x) &= \zeta(9)\left(\frac{340}{3} - \frac{146}{3}x - \frac{104}{9}x^2 - \frac{1}{4}x^3\right) + \zeta(3)^3\left(\frac{32}{3} + \frac{80}{3}x\right), \\ C_{10}(x) &= \zeta(7)\zeta(3)(x^3 + 39x^2 - 24x + 144) + \zeta(5)^2(x^3 - 6x^2 - 12x + 72), \\ C_{11}(x) &= \zeta(11)\left(-\frac{5}{32}x^4 - \frac{661}{60}x^3 - \frac{1623}{20}x^2 - 176x + 372\right) \\ &\quad + \zeta(5)\zeta(3)^2(80x^2 + 176x + 96) + \zeta^{\text{sv}}(3, 5, 3)\left(\frac{1}{5}x^3 + \frac{54}{5}x^2\right), \\ C_{12}(x) &= \zeta(9)\zeta(3)\left(\frac{5}{8}x^4 + \frac{107}{3}x^3 + 456x^2 - \frac{488}{3}x + \frac{1360}{3}\right) \\ &\quad + \zeta(7)\zeta(5)\left(\frac{11}{8}x^4 + \frac{47}{2}x^3 - 207x^2 - 216x + 432\right) + \zeta(3)^4(-16x^2 + \frac{272}{3}x + \frac{32}{3}), \\ C_{13}(x) &= \zeta(13)\left(-\frac{7}{64}x^5 - \frac{226501}{16800}x^4 - \frac{1283839}{8400}x^3 - \frac{1447393}{1400}x^2 - 618x + 1260\right) \\ &\quad + \zeta(7)\zeta(3)^2(x^4 + 34x^3 + 1236x^2 + 336x + 288) \end{aligned} $$
$$ 
\begin{aligned}
&+ \zeta(5)^2\zeta(3)\left(-\frac{31}{10}x^4 + \frac{256}{5}x^3 - \frac{1128}{5}x^2 + 336x + 288\right) \\
&+ \zeta^{\text{sv}}(5, 3, 5)\left(-\frac{157}{1400}x^4 - \frac{549}{350}x^3 - \frac{12339}{175}x^2\right) \\
&+ \zeta^{\text{sv}}(3, 7, 3)\left(\frac{5}{56}x^4 + \frac{59}{28}x^3 + \frac{747}{14}x^2\right).
\end{aligned}
$$

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 can be downloaded from [7]. 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, starting at the eleventh order. In the case of $C_{15}(x)$, the basis of single-valued multiple zeta values was found with the help of the program HYPERLOGPROCEDURES developed by Oliver Schnetz. In other cases, we have used the generators for the algebra of single-valued multiple zeta values given by Francis Brown in [8, §7.4]. 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 conditions in the companion paper [6]. Namely, we explicitly derived the results up to (and including) $C_{14}$ and proved 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.

## 2 The method of particular solutions

Because the eigenfunctions considered in this work are dihedral 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 was introduced by Conway [9] in 1961 and established by a famous paper of Fox, Henrici, and Moler [10] in 1966 on computing eigenvalues of L-shape domains.

The main idea of the MPS is to expand $\Psi$ for the spectral parameter $k = \sqrt{\lambda}$ as a Fourier-Bessel series (in polar coordinates)

$$ \Psi(k, r, \theta) = \sum_{v=1}^{\infty} c_v \psi_v(k, r, \theta), \qquad (4) $$

where

$$ \psi_v(k, r, \theta) = J_{mv}(k \cdot r) \cdot \begin{cases} \sin(m_v \theta) & \text{Odd eigenfunctions,} \\ \cos(m_v \theta) & \text{Even eigenfunctions,} \end{cases} \qquad (5) $$
<figure><img src="image_3.png" /><figcaption>Fig. 1 Fundamental region of regular polygons for fully symmetric eigenfunctions</figcaption></figure>

and $J_m$ denotes the Bessel function (see, for example, [11]). The basis functions $\psi_v$ have the useful property that they can serve as eigenfunctions along an unbounded wedge by a well-made choice of $m_v$. Consider for example a wedge with angle $\alpha$ with the vertex at the origin. Dirichlet boundary conditions impose that the function has to vanish along the boundary. If one chooses $m_v = v\pi/\alpha$ 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 $\psi_v$. To obtain the spectral parameter $k$ one truncates the expansion of $\Psi$ to some finite order $M$ and searches for the parameter $k$ for which the expansion vanishes on a discrete set of $N$ 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

$$ \begin{pmatrix} \psi_1(k, r_1, \theta_1) & \dots & \psi_M(k, r_1, \theta_1) \\ \vdots & \ddots & \vdots \\ \psi_1(k, r_N, \theta_N) & \dots & \psi_M(k, r_N, \theta_N) \end{pmatrix} \cdot \begin{pmatrix} c_1 \\ \vdots \\ c_M \end{pmatrix} = \begin{pmatrix} 0 \\ \vdots \\ 0 \end{pmatrix}, $$

which is square if the number of point-matching points is also chosen to be $M$. 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) [10]. By increasing $N$, one can obtain better approximations and hence better precision for the eigenvalue.

This approach has been applied by the second author [5, 12] who expanded from the vertex with angle $\beta/2$ (see Fig. 1) and applied point-matching to the edge $\partial\Omega_2$ to compute eigenvalues of polygons with $n \le 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 [13, 14] 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 [12].

We further remark that recent progress has been achieved in making the results of the MPS rigorous (i.e., with certified error bounds) [15, 16]. 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.

## 3 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 [5, 12] 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 the decrease of precision for polygons with more edges, we used the technique of domain decomposition. Domain decomposition is by no means new in the context of the MPS. The first application was done by J. Descloux and M. Tolley [17] in 1983 who decomposed domains into several subdomains $\Omega_k$ with eigenfunctions $f_k$ and minimized the quotient of the functionals

$$ \sum_{k<l}^{N} \int_{\Gamma_{kl}} (|f_k - f_l|^2 + |\nabla f_k - \nabla f_l|^2) \, ds, \quad (6) $$

and

$$ \sum_{k=1}^{N} \iint_{\Omega_k} |f_k|^2 \, dx \, dy, \quad (7) $$
where $\Gamma_{kl}$ denote straight intersecting boundaries. This method was further developed by Driscoll [18] in 1997 to compute eigenmodes of the famous isospectral drums of Gordon, Webb, and Wolpert [19]. Betcke [20] 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. (6), 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 [18]); however, we are unaware of any author proving the success of this method.

## 4 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

$$ \Psi_I^{[M_I]}(k, r, \theta) = \sum_{v=1}^{M_I} J_{m_v}(k \cdot r) \cdot \cos(m_v \theta), \quad (8) $$

and expands at the vertex with angle $\alpha/2$. To fulfill the Neumann boundary condition at angle $\theta = 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_v \cdot 0) = 0$ for all $m_v$. The second Neumann

<figure><img src="image_4.png" /><figcaption>Fig. 2 Domain decomposition of the fundamental region</figcaption></figure>
boundary condition is fulfilled if $\sin(m_v \frac{\alpha}{2}) = 0$ from which follows

$$m_v = \frac{2\pi}{\alpha}(\nu - 1). \qquad (9)$$

Similarly for the second domain one expands from the $\beta/2$-vertex (and places it at the origin). We denote this coordinate system with a tilde. The space of eigenfunctions now becomes

$$\Psi_{II}^{[M_{II}]}(k, \tilde{r}, \tilde{\theta}) = \sum_{\mu=1}^{M_{II}} J_{m_{\mu}}(k \cdot \tilde{r}) \cdot \sin(m_{\mu} \tilde{\theta}), \qquad (10)$$

with

$$m_{\mu} = \frac{\pi}{\beta}(2\mu - 1). \qquad (11)$$

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

$$\begin{aligned} \Psi_I(p_1) &= \Psi_{II}(\tilde{p}_1), \\ \vdots \\ \Psi_I(p_{N_p}) &= \Psi_{II}(\tilde{p}_{N_p}). \end{aligned}$$

2. The normal derivatives of both wave functions have to match along the intersecting line

$$\begin{aligned} \frac{\partial}{\partial x} \Psi_I(p_1) &= - \frac{\partial}{\partial \tilde{y}} \Psi_{II}(\tilde{p}_1) \\ \vdots \\ \frac{\partial}{\partial x} \Psi_I(p_{N_p}) &= - \frac{\partial}{\partial \tilde{y}} \Psi_{II}(\tilde{p}_{N_p}), \end{aligned}$$

where the minus sign appears because the expansions have opposite orientations along the intersection line. We computed the derivatives explicitly by differentiating the Bessel and sine functions.

3. The outer boundary condition of the second region has to be fulfilled

$$\begin{aligned} \frac{\partial}{\partial \tilde{x}} \Psi_{II}(\tilde{q}_1) &= 0, \\ \vdots \\ \frac{\partial}{\partial \tilde{x}} \Psi_{II}(\tilde{q}_{N_q}) &= 0. \end{aligned}$$

All these conditions can then be combined into a point-matching matrix
$$A = \begin{pmatrix}
\Psi_I(p_1) & -\Psi_{II}(\tilde{p}_1) \\
\vdots & \vdots \\
\Psi_I(p_{N_p}) & -\Psi_{II}(\tilde{p}_{N_p}) \\
\hline
\frac{\partial}{\partial x}\Psi_I(p_1) & \frac{\partial}{\partial y}\Psi_{II}(\tilde{p}_1) \\
\vdots & \vdots \\
\frac{\partial}{\partial x}\Psi_I(p_{N_p}) & \frac{\partial}{\partial y}\Psi_{II}(\tilde{p}_{N_p}) \\
\hline
0 & \frac{\partial}{\partial x}\Psi_{II}(\tilde{q}_1) \\
\vdots & \vdots \\
0 & \frac{\partial}{\partial x}\Psi_{II}(\tilde{q}_{N_p})
\end{pmatrix} \quad (12)$$

We chose the height of the intersecting line to be equal to the length of $\partial\Omega_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_v$). 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 $M_I = N/4$ and $M_{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.

## 5 Details on the numerical implementation

We implemented the algorithm in the JULIA programming language [21]. We made use of ARB [22] and its JULIA interface implemented in NEMO [23]. 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 [24] 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 [25]. 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
<figure><img src="image_5.png" /><figcaption>Fig. 3 Comparison of the eigenvalue precision for different polygons with 200 matching points</figcaption></figure>

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 computations 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 [26] which offers a 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 secant method) the result (and other significant parameters) were stored and the job re-queued. This approach allowed other users to occupy resources as well as 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 h 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 thread hours in total (the amount of physical core hours without hyper-threading is impossible to determine).

# 6 Derivation of the coefficients

To compute the coefficients of the 1/n-expansion, we used Richardson extrapolation as it is described in [27]. The computed eigenvalues $\lambda(n)$ for n-sided polygons can be

<figure><img src="image_6.png" /></figure>
expanded as a series

$$
\begin{align*}
\lambda(n) &= C_0 + C_1 \cdot n^{-1} + C_2 \cdot n^{-2} + \dots + C_N \cdot n^{-N} \\
\lambda(n+1) &= C_0 + C_1 \cdot (n+1)^{-1} + C_2 \cdot (n+1)^{-2} + \dots + C_N \cdot (n+1)^{-N} \\
\vdots \\
\lambda(n+N) &= C_0 + C_1 \cdot (n+N)^{-1} + C_2 \cdot (n+N)^{-2} + \dots + C_N \cdot (n+N)^{-N}
\end{align*}
$$

where *n* is the lowest computed polygon and *N* is the total number 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 numer-
ical estimates for the coefficients *C*ᵢ. Using this approach gives higher precision than
a standard polynomial fit as it was applied in [5]. 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 [28]) or by performing the computations with rational
arithmetic to not implement any rounding errors) but none of them yielded any better
results.

An empirical observation is that one seems to get the best results by using about 0.65*D* 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 [29]
to find closed-form guesses using the LLL algorithm [30]. For regular polygons with
Dirichlet boundary condition, we have also made use of the results of [6], which give
closed forms for the coefficients of x⁰ and x¹ (i.e., the parts of the coefficients Cᵢ
which are constant or linear in x), where x = j²₀,₁ as above, to cut down on the search
space.

**7 Convergence of the asymptotic series**

By approximating the growth rate of the coefficients in Eq. (3), one can attempt
to understand whether the asymptotic series converges. In general, there are three
possibilities:

1. The asymptotic series converges for all *n*;

2. The asymptotic series diverges for $n < n_{\text{crit}}$ and converges for $n \ge n_{\text{crit}}$ for some "critical" polygon with $n_{\text{crit}}$ edges;

3. The asymptotic series diverges for all *n*.

It was already observed by the second author in [5] that the expansion coefficients seem
to have at least exponential growth and alternating signs. He suggested the asymptotic
approximation formula

$$ \ln |C_i| \approx a \cdot i + b \cdot (1 - \exp(-c \cdot i)). \quad (13) $$

By fitting the computed data to this ansatz, it was thus conjectured in [5] that there exists some “critical” value $n_{crit}$ such that the series in Eq. (3) converges for $n \ge n_{crit}$ or diverges for $n < n_{crit}$. However, we found that the estimates for $n_{crit}$ continue to increase as more data is added to the fit. An alternative expansion that includes an $\ln-\ln$ term is given by

$$ \ln |C_i| \approx (a + \ln(\ln(i)))i + b(1 - \exp(-c \cdot i)), \quad (14) $$

with parameters

* $a = 2.00752 \pm 0.005682$

* $b = -254.721 \pm 2.557$

* $c = 7.70011e-2 \pm 7.46e-5$

results in a better fit with $\chi^2_{reduced} = 2.091$ (note that the error bounds are only related to the precision of the fit and do not provide real bounds). However, the parameter $a$ is still relatively dependent on the amount of data added to the fit ($a$ grows as more orders are included), making it impossible to discuss whether the asymptotic series converges or not. Table 1 shows the precision of the eigenvalue expansion formula, given in Eq. (3), for various $n$-gons at different truncation orders. Due to the increasing residuals at the cases $n = 3$ and $n = 4$, it becomes quite likely that these series diverge. Therefore, it seems unlikely that the asymptotic series converges for all $n$. However, it remains open whether the series is divergent for all $n$ or only for some $n < n_{crit}$.

<table><caption>**Table 1** Truncating Eq. (3) at some finite order $M$ gives eigenvalue estimates $\lambda(n)^{[M]}$</caption><thead><tr><th>n</th><th>M = 8</th><th>M = 10</th><th>M = 12</th><th>M = 14</th><th>M = 16</th></tr></thead><tbody><tr><td>3</td><td>0.0311437</td><td>0.053422</td><td>0.131248</td><td>0.270729</td><td>0.649124</td></tr><tr><td>4</td><td>4.2384e-04</td><td>3.43341e-03</td><td>3.82489e-03</td><td>4.29893e-03</td><td>5.77807e-03</td></tr><tr><td>5</td><td>2.17896e-04</td><td>3.52867e-04</td><td>2.37249e-04</td><td>1.69646e-04</td><td>1.46169e-04</td></tr><tr><td>6</td><td>6.16703e-05</td><td>5.30761e-05</td><td>2.4125e-05</td><td>1.19703e-05</td><td>7.17852e-06</td></tr><tr><td>7</td><td>1.88765e-05</td><td>1.05309e-05</td><td>3.46423e-06</td><td>1.2638e-06</td><td>5.58066e-07</td></tr><tr><td>8</td><td>6.47816e-06</td><td>2.57115e-06</td><td>6.41478e-07</td><td>1.79396e-07</td><td>6.07743e-08</td></tr><tr><td>9</td><td>2.46731e-06</td><td>7.37098e-07</td><td>1.4438e-07</td><td>3.19478e-08</td><td>8.56699e-09</td></tr><tr><td>10</td><td>1.02719e-06</td><td>2.40121e-07</td><td>3.79268e-08</td><td>6.807e-09</td><td>1.48094e-09</td></tr></tbody></table>

The entries of this table correspond to $|\lambda(n) - \lambda(n)^{[M]}|$. The cases $n = 3$ & $n = 4$ appear to be diverging
# 8 Further examples

## 8.1 Regular polygons with Neumann boundary condition

The presented algorithm works equally well for regular polygons with Neumann boundary conditions (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

$$ \lambda_1(n) \sim j_{1,1}^2 \sum_{k=0}^{\infty} \frac{C_k(j_{1,1}^2)}{n^k} \quad (15) $$

for the Neumann case with coefficients (here $x = j_{1,1}^2$)

$$ C_0(x) = 1, $$

$$ C_1(x) = 0, $$

$$ C_2(x) = 0, $$

$$ C_3(x) = 0, $$

$$ C_4(x) = 0, $$

$$ C_5(x) = -4\zeta(5)x, $$

$$ C_6(x) = 0, $$

$$ C_7(x) = \zeta(7)(-2x^2 - 22x), $$

$$ C_8(x) = 0, $$

$$ C_9(x) = \zeta(9)(-\frac{3}{2}x^3 - \frac{52}{3}x^2 - \frac{296}{3}x), $$

$$ C_{10}(x) = \zeta(5)^2(2x^3 + 16x^2), $$

$$ C_{11}(x) = \zeta(11)(-\frac{5}{4}x^4 - \frac{215}{12}x^3 - 108x^2 - 418x), $$

$$ C_{12}(x) = \zeta(7)\zeta(5)(4x^4 + 46x^3 + 176x^2), $$

$$ C_{13}(x) = -\zeta(13)\left(\frac{35}{32}x^5 + \frac{4043}{150}x^4 + \frac{1699}{12}x^3 + 590x^2 + 1732x\right) - \frac{48}{5}\zeta(5)^2\zeta(3)x^4 - \frac{2}{25}\zeta^{\text{sv}}(5, 3, 5)x^4, $$

$$ C_{14}(x) = \zeta(9)\zeta(5)\left(\frac{15}{4}x^5 + \frac{826}{9}x^4 + \frac{920}{3}x^3 + \frac{2368}{3}x^2\right) + \zeta(7)^2\left(\frac{9}{4}x^5 + 18x^4 + \frac{385}{2}x^3 + 484x^2\right), $$

$$ C_{15}(x) = -\zeta(15)\left(\frac{63}{64}x^6 - \frac{543919}{6720}x^5 + \frac{16651961}{18900}x^4 + \frac{33557}{36}x^3 + \frac{15014}{5}x^2 + \frac{35494}{5}x\right) - \zeta(5)^3\left(\frac{20}{3}x^5 + \frac{1288}{15}x^4 + 64x^3\right) - \zeta^{\text{sv}}(5, 3, 7) + 336\zeta(7)\zeta(5)\zeta(3)\left(\frac{1}{14}x^5 - \frac{13}{35}x^4\right), $$

$$ C_{16}(x) = \zeta(11)\zeta(5)\left(\frac{7}{2}x^6 + \frac{16993}{120}x^5 + \frac{7088}{5}x^4 + 1714x^3 + 3344x^2\right) + \zeta(9)\zeta(7)\left(\frac{9}{2}x^6 + \frac{6409}{72}x^5 + \frac{1231}{9}x^4 + \frac{6836}{3}x^3 + \frac{13024}{3}x^2\right) + \zeta(5)^2\zeta(3)^2(12x^5 - 192x^4) + \zeta^{\text{sv}}(3, 5, 3)\zeta(5)\left(\frac{6}{5}x^5 - \frac{96}{5}x^4\right). $$
Our results thus indicate that the appearance of single-valued MZVs is not restricted to the Dirichlet eigenvalues (although the results for the Neumann boundary condition have not been formally proven yet).

## 8.2 Regular star polygons

In order to study whether the expandability of eigenvalues also applies to 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 Schlafi 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 ($\pi/n$, $\pi(n-2q)/2n$, $\pi(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 [12] 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_1 \approx N_{II} \approx N_{III} \approx 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 \le n \le 263$ to almost 300 digits precision and in the cases $3 \le q \le 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

$$ 
\begin{aligned}
C_0 &= 1, \\
C_1 &= 0, \\
C_2 &= 0, \\
C_3 &= \zeta(3)(14q(q-1)+4),
\end{aligned}
$$

<figure><img src="image_7.png" /><figcaption>Fig. 4 Regular star polygons with 11 edges and different densities. The figure was created using a MATH-EMATICA package by Michael Schreiber [31]</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 5 Domain decomposition of the fundamental region of regular star polygons</figcaption></figure>

$$C_4 = 0,$$
$$C_5 = \zeta(5)(62q(q-1)(q^2-q+1) + 12 - x \cdot \frac{31}{4}(q^2-q+\frac{8}{31})),$$

$$C_6 = \zeta(3)^2(2(7q(q-1)+2)^2 + x \cdot \frac{63}{4}(q^2-q+\frac{16}{63})),$$

$$C_7 = \zeta(7)(254q(q-1)(q^2-q+1)^2 + 36) - x \cdot A_7q(q-1)(2q-1)^2 \\ + x \cdot \zeta(7)(\frac{809}{12}q(q-1)(q^2-q+1) - \frac{1571}{16}q(q-1) - 12) \\ - x^2 \cdot \zeta(7)(\frac{127}{64}q(q-1) + \frac{1}{2}),$$

$$C_8 = \zeta(3)\zeta(5)(28q(q-1)+8)(31q(q-1)(q^2-q+1)+6) \\ + x \cdot \zeta(3)\zeta(5)(110q(q-1)(q^2-q+\frac{95}{176})+8) \\ - x \cdot 28\zeta(3)(\zeta(2)\zeta(3) - 2\mathrm{Li}_{3,1}(1,1,1))q(q-1)(2q-1)^2 \\ + x^2 \cdot \zeta(3)\zeta(5)(\frac{255}{32}q(q-1)+2).$$

Here the number $A_7$ is defined as

$$A_7 = \frac{56}{9}\mathrm{Li}_{2,2,3}(-1, -1, -1) + 8\mathrm{Li}_{2,2,3}(-1, 1, -1) + \frac{7435}{288}\zeta(5, 2) + \frac{775}{144}\zeta(4, 3).$$

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.

## 8.3 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 a radius

$$r(\theta) = R + d \cdot \cos(n \cdot \theta) \quad (16)$$

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
<figure><img src="image_9.png" /><figcaption>Fig. 6 Smooth star shape with $n = 5$ and $d = 0.3$</figcaption></figure>

been investigated by Laura in 1964 [32] and Wagner in 1971 [33] who derived the approximate formula

$$ \lambda(n, d) = c^2 \cdot \left[ 1 + 0.125 c \cdot \left( \frac{J_2(c)}{J_1(c)} + \frac{2(J_{n-1}(c) - 2J_{n+1}(c))}{J_n(c)} \right) \cdot d^2 \right], \quad (17) $$

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, [34, 35]. In our case, we investigate the behavior of smooth star shape domains with $d = (1/n)^m$ where $m \in \mathbb{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 Fig. 7. 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_{0,1}^2$. For example for $m = 2$ one obtains the rational polynomials

$$ C_0(x) = 1, $$

$$ C_1(x) = 0, $$

$$ C_2(x) = 0, $$
<figure><img src="image_10.png" /><figcaption>Fig. 7 Hypocycloid with $n = 10$</figcaption></figure>

$$C_3(x) = 1,$$

$$C_4(x) = 1/2,$$

$$C_5(x) = -(2x + 1)/4,$$

$$C_6(x) = (2x + 3)/4,$$

$$C_7(x) = -(12x^2 + 18x - 119)/96,$$

$$C_8(x) = (48x^2 - 75x - 8)/96,$$

$$C_9(x) = -(144x^3 + 3393x^2 - 84x + 167)/2304,$$

$$C_{10}(x) = (384x^3 + 2499x^2 + 712x + 1203)/768,$$

$$C_{11}(x) = -(14400x^4 + 926100x^3 + 2283120x^2 + 525440x - 390979)/368640,$$

$$C_{12}(x) = (46080x^4 + 933075x^3 + 1041030x^2 - 95940x - 63551)/92160.$$

(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.

## 8.4 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 [36]. Hypocycloids can be parameterized by

$$ 
\begin{aligned}
x(\theta) &= r \cdot (n-1) \cdot \cos(\theta) + r \cdot \cos((n-1) \cdot \theta) \\
y(\theta) &= r \cdot (n-1) \cdot \sin(\theta) - r \cdot \sin((n-1) \cdot \theta)
\end{aligned}
$$

where $r = R/n$ is the radius of the smaller (rotating) circle. In the limit as $n \to \infty$ 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 cusp) is therefore the general Fourier-Bessel expansion with integer coefficients $m_v$

$$ \Psi_{II}^{[M]}(k, r, \theta) = \sum_{m_v=1}^{M_{II}} J_{m_v}(k \cdot r) \cdot \cos(m_v \cdot \theta). \quad (18) $$

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 computed each hypocycloid eigenvalue from 5 to 100 edges for 1 week on a single thread, which resulted in around 65 digits of precision per eigenvalue. If one fits the eigenvalues on a 1/n-polynomial-expansion, the expansion coefficients are converging toward constant values which are given by

$$ 
\begin{aligned}
C_0 &= 1 \\
C_1 &= 0.372285278452919538950832494889900 \\
C_2 &= 1.16087789045807729220949677611 \\
C_3 &= 1.8723005009420124517917316
\end{aligned}
$$

<figure><img src="image_11.png" /><figcaption>Fig. 8 Decomposition of the fundamental region of a hypocycloid</figcaption></figure>
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 [37]), we found indications that the expansion coefficients $C_i$ are polynomials of degree $\lfloor \frac{i-1}{2} \rfloor$ in the limiting eigenvalue $x_m = j_{0,m}^2$. Namely, by computing the first four fully symmetric eigenvalues of several hypocycloids, we found evidence that

$$C_0 = 1 \cdot (x_m)^0$$

$$C_1 = 0.372285278452919538950832494889900 \cdot (x_m)^0$$

$$C_2 = 1.16087789045807729220949677611 \cdot (x_m)^0$$

$$C_3 = 2.4397482604 \cdot (x_m)^0 - 0.0981202685 \cdot (x_m)^1$$

$$C_4 = 4.9228013 \cdot (x_m)^0 - 0.088251328 \cdot (x_m)^1$$

However, it remains unclear whether these numbers can be expressed in closed form.

## 9 Conclusion

We have demonstrated how the method of particular solutions combined with domain decomposition can be efficiently applied to compute Laplacian eigenvalues of several two-dimensional shapes with dihedral symmetry to very high precision. Using this numerical data, we have computed two additional eigenvalue expansion coefficients of regular polygons with Dirichlet boundary condition as well as the first 17 coefficients for regular polygons with Neumann boundary condition and showed that single-valued multiple zeta values appear in the expansions of both of these examples. Moreover, we have also shown that Riemann zeta values and alternating MZVs appear in the eigenvalue expansion coefficients of more general regular star shapes, together with a formula relating the coefficients of star shapes with different densities. Additionally, we constructed examples of star shapes with sinusoidal boundaries and found that the expansion coefficients are defined over $\mathbb{Q}[x]$ which indicates that the $1/n$ expandability of eigenvalues is not restricted to polygonal shapes only.

## Appendix 1. Brief introduction to multiple zeta values

The Riemann zeta function

$$\zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s} \qquad (A1)$$

is known to converge for all $s \in \mathbb{C}$ that satisfy $\text{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 $\zeta(s)$ at natural numbers $s \in \mathbb{N}$, $s > 1$, which are sometimes simply called zeta values. For even zeta values, Euler has famously shown that

$$ \zeta(2n) = \frac{(-1)^{n+1} B_{2n} (2\pi)^{2n}}{2(2n)!}, \qquad (A2) $$

where $B_n$ are Bernoulli numbers. For odd zeta values much less is known. It has been proven by Apéry in 1979 that $\zeta(3)$ is irrational. His proof has not yet been extended to any other odd zeta values. It is known that at least one of $\zeta(5)$, $\zeta(7)$, $\zeta(9)$, $\zeta(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 $\pi$, $\zeta(3)$, $\zeta(5)$, ..., $\zeta(2n + 1)$, ... have any algebraic relations among them (altough it is widely believed that there are none). As we have seen in Sect. 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 $\zeta(7)\zeta(3)$ and $\zeta(5)^2$, and we have $7 + 3 = 2 \cdot 5 = 10$. This number is called the weight of a zeta product. Note that the product of two zeta values can be written as

$$ \zeta(a)\zeta(b) = \sum_{n=1}^{\infty} \frac{1}{n^a} \cdot \sum_{m=1}^{\infty} \frac{1}{m^b} = \sum_{n,m=1}^{\infty} \frac{1}{n^a m^b}. \qquad (A3) $$

We can split the last sum into three terms

$$ \sum_{n,m=1}^{\infty} \frac{1}{n^a m^b} = \left( \sum_{0<n<m} + \sum_{0<n=m} + \sum_{0<m<n} \right) \frac{1}{n^a m^b}. \qquad (A4) $$

One of the terms is given as a zeta value

$$ \sum_{n=m=1}^{\infty} \frac{1}{n^a m^b} = \sum_{n=1}^{\infty} \frac{1}{n^{a+b}} = \zeta(a+b). \qquad (A5) $$

The remaining “crossing” terms are the so-called multiple zeta values $\zeta(a, b)$ and $\zeta(b, a)$. More generally, the multiple zeta values (MZVs) are defined as

$$ \zeta(s_1, \dots, s_k) = \sum_{0 < n_1 < \dots < n_k} \frac{1}{n_1^{s_1} \dots n_k^{s_k}}. \qquad (A6) $$

The decomposition

$$ \zeta(a)\zeta(b) = \zeta(a, b) + \zeta(a+b) + \zeta(b, a) \qquad (A7) $$

is sometimes called the Nielsen reflection formula. Multiple zeta values satisfy a wide variety of interesting relations. For instance, the $\mathbb{Q}$-linear span of MZVs inside $\mathbb{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

$$ \zeta^{\text{sv}}(2n + 1) = 2\zeta(2n + 1). \qquad (A8) $$

We will not cover the theory of single-valued MZVs in this brief introduction but instead, only give an expression for the single-valued MZVs that appear in our expansion formulas

$$ \zeta^{\text{sv}}(3, 5, 3) = 2\zeta(3, 5, 3) - 2\zeta(3)\zeta(5, 3) - 10\zeta(3)^2\zeta(5), $$

$$ \zeta^{\text{sv}}(5, 3, 5) = 2\zeta(5, 3, 5) - 22\zeta(5)\zeta(5, 3) - 120\zeta(5)^2\zeta(3) - 10\zeta(5)\zeta(8), $$

$$ \zeta^{\text{sv}}(3, 7, 3) = 2\zeta(3, 7, 3) - 2\zeta(3)\zeta(7, 3) - 28\zeta(3)^2\zeta(7) - 24\zeta(5)\zeta(5, 3) \\ - 144\zeta(5)^2\zeta(3) - 12\zeta(5)\zeta(8), $$

$$ \zeta^{\text{sv}}(5, 3, 7) = 2\zeta(5, 3, 7) - 28\zeta(3, 5)\zeta(7) - 12\zeta(3, 7)\zeta(5) - 336\zeta(3)\zeta(5)\zeta(7) \\ - 78\zeta(5)^3 - 6\zeta(10)\zeta(5) - 28\zeta(8)\zeta(7) - 34\zeta(6)\zeta(9) + 110\zeta(4)\zeta(11) \\ + 1001\zeta(2)\zeta(13), $$

$$ \zeta^{\text{sv}}(3, 3, 9) = 2\zeta(3, 3, 9) + 12\zeta(3, 7)\zeta(5) + 30\zeta(3, 5)\zeta(7) + 272\zeta(5)^3 - 27\zeta(3)^2\zeta(9) \\ + 1209\zeta(13)\zeta(2) + 252\zeta(11)\zeta(4) + 29\zeta(9)\zeta(6) + 6\zeta(7)\zeta(8) \\ + 318\zeta(3)\zeta(5)\zeta(7), $$

$$ \zeta^{\text{sv}}(1, 1, 3, 4, 6) = 2\zeta(1, 1, 3, 4, 6) + 2\zeta(1, 1, 4, 6)\zeta(3) - 12\zeta(3)\zeta(3, 7)\zeta(2) \\ + 4\zeta(3, 3, 7)\zeta(2) + 20\zeta(3, 3, 5)\zeta(4) - \frac{13133}{56}\zeta(5)^3 - \frac{2}{3}\zeta(3)^5 \\ - \frac{23}{2}\zeta(3)^2\zeta(5)\zeta(4) + 21\zeta(3, 5)\zeta(5)\zeta(2) - 16\zeta(3)\zeta(5)^2\zeta(2) \\ - \frac{481}{10}\zeta(3, 5)\zeta(7) + \frac{349}{45}\zeta(7)\zeta(8) - \frac{118175}{336}\zeta(5)\zeta(10) - \frac{2927}{9}\zeta(3)^2\zeta(9) \\ - \frac{56717}{60}\zeta(13)\zeta(2) + \frac{84185}{72}\zeta(9)\zeta(6) - \frac{27199}{12}\zeta(11)\zeta(4) - \frac{28}{5}\zeta(2)\zeta(3, 5, 5) \\ + \frac{58}{9}\zeta(3)\zeta(3, 9) + \frac{3992923}{49752}\zeta(3)\zeta(12) + 3\zeta(3)^3\zeta(6) - \frac{145}{56}\zeta(5)\zeta(3, 7) \\ - 8\zeta(3)\zeta(3, 5)\zeta(4) - \frac{12509}{24}\zeta(3)\zeta(5)\zeta(7) - 52\zeta(3)^2\zeta(7)\zeta(2). $$

Some of the expressions, we give also involve alternating multiple zeta values which are special values of multiple polylogarithms

$$ \mathrm{Li}_{m_1, \dots, m_k}(z_1, \dots, z_k) = \sum_{0 < n_1 < \dots < n_k} \frac{z_1^{n_1} \dots z_k^{n_k}}{n_1^{m_1} \dots n_k^{m_k}}, \qquad (A9) $$

when all $z_i$ are set to $\pm 1$.
Remaining coefficients for regular polygons with Dirichlet boundary condition

$$
\begin{aligned}
C_{14}(x) ={}& \zeta(11)\zeta(3)(1488 - 664x + \frac{32902}{5}x^2 + \frac{20381}{30}x^3 + \frac{10169}{240}x^4 + \frac{7}{16}x^5) \\
& + \zeta(9)\zeta(5)(1360 - \frac{2824}{3}x - \frac{7306}{3}x^2 + \frac{1300}{9}x^3 + \frac{467}{8}x^4 + x^5) \\
& + \zeta(7)^2(648 - 540x - \frac{3627}{2}x^2 + \frac{483}{4}x^3 + \frac{175}{16}x^4 + \frac{9}{16}x^5) \\
& + \zeta(5)\zeta(3)^3(128 + \frac{2752}{3}x + \frac{3664}{3}x^2 + 40x^3) \\
& + \zeta^{\text{sv}}(3, 5, 3)\zeta(3)x^2 \frac{1296+28x-x^2}{5},
\end{aligned}
$$

$$
\begin{aligned}
C_{15}(x) = & \left(-\frac{96}{11}x^3 - \frac{864}{111}x^2\right)\zeta^{\text{sv}}(1, 1, 3, 4, 6) \\
& + \left(-\frac{7}{72}x^5 - \frac{71}{36}x^4 - \frac{422}{111}x^3 - \frac{6372}{111}x^2\right)\zeta^{\text{sv}}(3, 3, 9) \\
& + \left(-\frac{19}{168}x^5 - \frac{2173}{1680}x^4 - \frac{13452}{385}x^3 - \frac{282042}{385}x^2\right)\zeta^{\text{sv}}(5, 3, 7) \\
& + \left(-\frac{21}{256}x^6 + \frac{2008799}{5040}x^5 + \frac{644330219}{100800}x^4 + \frac{12325006817}{103950}x^3\right) \\
& + \left(\frac{26712690853}{11550}x^2 - \frac{10852}{5}x + \frac{21844}{5}\right)\zeta(15) \\
& + \left(-\frac{15}{4}x^5 - \frac{2183}{18}x^4 - 4842x^3 - \frac{158368}{3}x^2 + \frac{2128}{3}x + \frac{2720}{3}\right)\zeta(9)\zeta(3)^2 \\
& + \left(-\frac{133}{12}x^5 + \frac{646}{15}x^4 - \frac{274916}{55}x^3 - \frac{4524336}{55}x^2 + 864x + 1728\right)\zeta(7)\zeta(5)\zeta(3) \\
& + \left(-3x^5 - \frac{511}{10}x^4 - \frac{420268}{165}x^3 - \frac{2076696}{55}x^2 + 144x + 288\right)\zeta(5)^3 \\
& + \left(-\frac{1664}{55}x^3 - \frac{62528}{165}x^2 + \frac{3136}{15}x + \frac{128}{15}\right)\zeta(3)^5,
\end{aligned}
$$

$$
\begin{aligned}
C_{16}(x) = & \left(\frac{21}{64}\zeta(13)\zeta(3) + \frac{49}{64}\zeta(11)\zeta(5) + \frac{29}{32}\zeta(9)\zeta(7)\right)x^6 \\
& + \left(\frac{446687}{8400}\zeta(13)\zeta(3) + \frac{42067}{480}\zeta(11)\zeta(5) + \frac{15073}{288}\zeta(9)\zeta(7) - 2\zeta(7)\zeta(3)^3\right) \\
& + \frac{32}{5}\zeta(5)^2\zeta(3)^2 + \left(\frac{67}{350}\zeta^{\text{sv}}(5, 3, 5) - \frac{1}{7}\zeta^{\text{sv}}(3, 7, 3)\right)\zeta(3) - \frac{1}{20}\zeta^{\text{sv}}(3, 5, 3)\zeta(5))x^5 \\
& + \left(\frac{3100857}{2800}\zeta(13)\zeta(3) + \frac{139433}{80}\zeta(11)\zeta(5) + \frac{11185}{72}\zeta(9)\zeta(7) - 20\zeta(7)\zeta(3)^3\right) \\
& - \frac{2584}{5}\zeta(5)^2\zeta(3)^2 - \left(\frac{1139}{350}\zeta^{\text{sv}}(5, 3, 5) + \frac{51}{28}\zeta^{\text{sv}}(3, 7, 3)\right)\zeta(3) - \frac{47}{10}\zeta^{\text{sv}}(3, 5, 3)\zeta(5))x^4 \\
& + \left(\frac{559169}{60}\zeta(13)\zeta(3) + \frac{2224}{5}\zeta(11)\zeta(5) + \frac{3721}{6}\zeta(9)\zeta(7) + 988\zeta(7)\zeta(3)^3\right) \\
& - 2808\zeta(5)^2\zeta(3)^2 - \left(\frac{522}{5}\zeta^{\text{sv}}(5, 3, 5) - 59\zeta^{\text{sv}}(3, 7, 3)\right)\zeta(3) + \frac{552}{5}\zeta^{\text{sv}}(3, 5, 3)\zeta(5))x^3 \\
& + \left(\frac{16964568}{175}\zeta(3)\zeta(13) - 26454\zeta(5)\zeta(11) - 40787\zeta(9)\zeta(7) + 18056\zeta(7)\zeta(3)^3\right) \\
& + \left(\frac{150864}{5}\zeta(5)^2\zeta(3)^2 - \left(\frac{214488}{175}\zeta^{\text{sv}}(5, 3, 5) - \frac{8964}{7}\zeta^{\text{sv}}(3, 7, 3)\right)\zeta(3)\right)x^2 \\
& + (-2424\zeta(13)\zeta(3) - 3480\zeta(11)\zeta(5) - 4184\zeta(9)\zeta(7) + 2240\zeta(7)\zeta(3)^3 + 3360\zeta(5)^2\zeta(3)^2)x \\
& + (5040\zeta(13)\zeta(3) + 4464\zeta(11)\zeta(5) + 4080\zeta(9)\zeta(7) + 384\zeta(7)\zeta(3)^3 + 576\zeta(5)^2\zeta(3)^2).
\end{aligned}
$$

Remaining coefficients for smooth star shapes with m = 2

$$
C_{13}(x) = -(1209600x^5 + 169412175x^4 + 1584305000x^3 + 991331400x^2 \\
- 67630780x - 22055334)/44236800,
$$

$$
C_{14}(x) = (13271040x^5 + 610174215x^4 + 3119672640x^3 + 1238016540x^2 \\
+ 7880224x + 62326230)/26542080,
$$

$$
C_{15}(x) = -(457228800x^6 + 118856311875x^5 + 2664517542050x^4 + 8258849383400x^3 \\
+ 2125616765440x^2 + 77434447468x - 10887717978)/22295347200,
$$

$$
C_{16}(x) = (33443020800x^6 + 2924868777075x^5 + 38003282808750x^4 + 76583141218500x^3 \\
+ 12609060010540x^2 + 18024777715x - 83085435396)/66886041600,
$$

<figure><img src="image_12.png" /></figure>
$$
\begin{aligned}
C_{17}(x) = & -(1810626048000x^7 + 785618449756875x^6 + 34479052407093000x^5 \\
& + 286806810669211000x^4 + 393746334459690800x^3 + 41833381428570655x^2 \\
& - 382662625041172x - 164371240997550)/112368549888000,
\end{aligned}
$$

$$
\begin{aligned}
C_{18}(x) = & (120394874880000x^7 + 17871670961353125x^6 + 468925800452550000x^5 \\
& + 2662356634044423750x^4 + 2565871095884593000x^3 + 179378699791862375x^2 \\
& - 584375985329753x + 732859733409750)/240789749760000,
\end{aligned}
$$

$$
\begin{aligned}
C_{19}(x) = & -(21184324761600000x^8 + 14232900142095075000x^7 + 1079843748016254200000x^6 \\
& + 18704185436989446975000x^5 + 75684321204286594467500x^4 \\
& + 52194215245123302758750x^3 + 2424740991307168439892x^2 + 8697289229328086768x \\
& + 806072685794305875)/1618107118387200000,
\end{aligned}
$$

$$
\begin{aligned}
C_{20}(x) = & (866843099136000000x^8 + 201696166234942312500x^7 + 9330236455363230750000x^6 \\
& + 113623492626641309812500x^5 + 337606293009361365696875x^4 \\
& + 168745237201472986724250x^3 + 5201634554633731393011x^2 \\
& + 7062388870424537486x - 2948155439034622500)/1733686198272000000,
\end{aligned}
$$

$$
\begin{aligned}
C_{21}(x) = & -(1779483279974400000000x^9 + 1750724028590009405625000x^8 \\
& + 210914990771238755890000000x^7 + 6570081043597798567754375000x^6 \\
& + 586977123215363713009888900000x^5 + 130785026894321633100311176125x^4 \\
& + 47808464146250191120470311872x^3 + 976817476078605785634707984x^2 \\
& - 709576341339684500163424x - 466588962879561812430000)/163105197533429760000000,
\end{aligned}
$$

$$
\begin{aligned}
C_{22}(x) = & (214075571762626560000000x^9 + 736244702627870619077343750x^8 \\
& + 54826137837085868012812500000x^7 + 1224633170486034514394552343750x^6 \\
& + 8269403144636044633016590206250x^5 + 14026490060106825798594234715375x^4 \\
& + 3774510137070113767485677126049x^3 + 51231956845133382896804013256x^2 \\
& - 34284798535733782456620476x \\
& + 15568092692374386637462500)/4281511435252531200000000,
\end{aligned}
$$

$$
\begin{aligned}
C_{23}(x) = & -(586994590598595379200000000x^{10} + 8099489427400822789150312500000x^9 \\
& + 1456333964071948513322434125000000x^8 + 74056060297895723139610998046875000x^7 \\
& + 1237529454152663240037004207720702500x^6 + 6452892881074851056674758533001245925x^5 \\
& + 8424983842296738700856453921666897340x^4 + 1676926979683976178266295659995421360x^3 \\
& + 15156930607344875353880119232561056x^2 + 3433373767917417187848282765632x \\
& + 1242040349979959360991463800000)/632978650587734212608000000000.
\end{aligned}
$$

Coefficients for smooth star shapes with m = 3

$C_0(x) = 1,$

$C_1(x) = 0,$

$C_2(x) = 0,$

$C_3(x) = 0,$

$C_4(x) = 0,$

$C_5(x) = 1,$

$C_6(x) = 1/2,$
$$C_7(x) = -x/2,$$

$$C_8(x) = x/2,$$

$$C_9(x) = -(x^2 + 4x + 2)/8,$$

$$C_{10}(x) = (2x^2 + 2x + 3)/4,$$

$$C_{11}(x) = -(x + 1)(x^2 + 21x - 18)/16,$$

$$C_{12}(x) = (16x^3 + 104x^2 - 25x + 10)/32,$$

$$C_{13}(x) = -(30x^4 + 1920x^3 + 5547x^2 - 192x - 88)/768,$$

$$C_{14}(x) = (384x^4 + 7776x^3 + 11523x^2 - 192x - 304)/768,$$

$$C_{15}(x) = -(168x^5 + 23520x^4 + 224715x^3 + 184920x^2 - 224x + 64)/6144,$$

$$C_{16}(x) = (3072x^5 + 141312x^4 + 758013x^3 + 367512x^2 + 5696x + 8192)/6144.$$

$$C_{17}(x) = (-12096x^6 - 3144960x^5 - 70977933x^4 - 235020744x^3 - 69971328x^2 - 927744x + 602368)/589824,$$

$$C_{18}(x) = (98304x^6 + 8601600x^5 + 113600385x^4 + 246217560x^3 + 46297920x^2 - 15360x + 88832)/196608.$$

Coefficients for smooth star shapes with $m = 4$

$$C_0(x) = 1,$$

$$C_1(x) = 0,$$

$$C_2(x) = 0,$$

$$C_3(x) = 0,$$

$$C_4(x) = 0,$$

$$C_5(x) = 0,$$

$$C_6(x) = 0,$$

$$C_7(x) = 1,$$

$$C_8(x) = 1/2,$$

$$C_9(x) = -x/2,$$

$$C_{10}(x) = x/2,$$

$$C_{11}(x) = -x(x + 4)/8,$$

$$C_{12}(x) = x(x + 1)/2,$$

$$C_{13}(x) = -(x^3 + 22x^2 + 8x + 4)/16,$$

$$C_{14}(x) = (2x^3 + 13x^2 + 2x + 3)/4,$$

$$C_{15}(x) = -(5x^4 + 320x^3 + 912x^2 + 24x - 144)/128,$$

$$C_{16}(x) = (16x^4 + 324x^3 + 480x^2 - 25x + 10)/32,$$

$$C_{17}(x) = -x(7x^4 + 980x^3 + 9360x^2 + 7929x - 64)/256,$$
$$C_{18}(x) = x(128x^4 + 5888x^3 + 31584x^2 + 16065x - 64)/256,$$

$$C_{19}(x) = (-126x^6 - 32760x^5 - 739344x^4 - 2453067x^3 - 773208x^2 + 1536x + 704)/6144,$$

$$C_{20}(x) = (3072x^6 + 268800x^5 + 3550080x^4 + 7730685x^3 + 1545624x^2 - 1536x - 2432)/6144.$$

**Acknowledgements** D.B. would like to thank Oliver Freyermuth for his help in carrying out the computations on a computing cluster, Fredrik Johansson for assistance with Arb, and Plamen Koev for discussions on possible improvements of the approach of Sect. 6. D.R. would like to thank Erik Panzer for his help with finding a good basis for single-valued MZVs in weight 16. D.R. would also like to thank Steven Charlton for his great help with calculations involving alternating MZVs. The authors would also like to thank the anonymous referees for their valuable comments. D.B. acknowledges financial support from the Bonn-Cologne Graduate School of Physics and Astronomy honors branch.

**Funding** Open Access funding enabled and organized by Projekt DEAL.

## Declarations

**Conflict of interest** The authors declare no competing interests.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

1. Pólya, G., Szegő, G.: Isoperimetric inequalities in mathematical physics. Ann. Math. Stud. (1951)

2. Molinari, L.: On the ground state of regular polygonal billiards. J. Phys. A Math. Gen. 30(18), 6517–6524 (1997). https://doi.org/10.1088/0305-4470/30/18/025

3. Grinfeld, P., Strang, G.: Laplace eigenvalues on regular polygons: a series in 1/n. J. Math. Anal. Appl. 385(1), 135–149 (2012). https://doi.org/10.1016/j.jmaa.2011.06.035

4. Boady, M.: Applications of symbolic computation to the calculus of moving surfaces. Phd thesis, Drexel University (2016)

5. Jones, R.S.: The fundamental Laplacian eigenvalue of the regular polygon with Dirichlet boundary conditions (2017). arXiv:1712.06082

6. Berghaus, D., Georgiev, B., Monien, H., Radchenko, D.: On Dirichlet eigenvalues of regular polygons (2021). arXiv:2103.01057

7. Berghaus, D., Jones, R.S., Monien, H., Radchenko, D.: Eigenvalue data (2022). https://github.com/David-Berghaus/master-thesis-data. Accessed 23 Sept 2022
8. Brown, F.: Single-valued motivic periods and multiple zeta values. Forum Math. Sigma. **2**, 25 (2014). https://doi.org/10.1017/fms.2014.18

9. Conway, H.D.: The bending, buckling, and flexural vibration of simply supported polygonal plates by point-matching. J. Appl. Mech. **28**(2), 288–291 (1961). https://doi.org/10.1115/1.3641670. https://asmedigitalcollection.asme.org/appliedmechanics/articlepdf/28/2/288/5444109/2881.pdf

10. Fox, L., Henrici, P., Moler, C.: Approximations and bounds for eigenvalues of elliptic operators. SIAM J. Numer. Anal. (1967). https://doi.org/10.1137/0704008

11. Vekua, I.N.: New methods for solving elliptic equations. North-Holland Series in Applied Mathematics and Mechanics, vol. 1, p. 358. North-Holland Publishing Co., Amsterdam; Interscience Publishers John Wiley & Sons, Inc., New York, ??? (1967)

12. Jones, R.S.: Computing ultra-precise eigenvalues of the Laplacian within polygons. Adv. Comput. Math. (2017). https://doi.org/10.1007/s10444-017-9527-y

13. Betcke, T., Trefethen, L.N.: Reviving the method of particular solutions. SIAM Review. **47**(3), 469–491 (2005)

14. Betcke, T.: The generalized singular value decomposition and the method of particular solutions. SIAM J. Sci. Comput. **30**(3), 1278–1295 (2008). https://doi.org/10.1137/060651057

15. Gómez-Serrano, J., Orriols, G.: Any three eigenvalues do not determine a triangle (2019). https://doi.org/10.48550/ARXIV.1911.06758. arXiv:1911.06758

16. Dahne, J., Salvy, B.: Computation of tight enclosures for Laplacian eigenvalues. SIAM J. Sci. Comput. **42**(5), 3210–3232 (2020). https://doi.org/10.1137/20M1326520

17. Descloux, J., Tolley, M.C.: An accurate algorithm for computing the eigenvalues of a polygonal membrane. Comput. Methods Appl. Mech. Eng. **39**(1), 37–53 (1983)

18. Driscoll, T.A.: Eigenmodes of isospectral drums. SIAM Rev., 39(1), 1–17 (1997). https://doi.org/10.1137/S0036144595285069

19. Gordon, C., Webb, D.L., Wolpert, S.: One cannot hear the shape of a drum. Bull. Amer. Math. Soc. **27**, 134–138 (1992). arXiv:math/9207215

20. Betcke, T.: A GSVD formulation of a domain decomposition method for planar eigenvalue problems. IMA Journal of Numerical Analysis. **27**(3), 451–478 (2006). https://doi.org/10.1093/imanum/drl030. https://academic.oup.com/imajna/article-pdf/27/3/451/2230584/drl030.pdf

21. Bezanson, J., Karpinski, S., Shah, V.B., Edelman, A.: Julia: A fast dynamic language for technical computing (2012). arXiv:1209.5145

22. Johansson, F.: Arb: efficient arbitrary-precision midpoint-radius interval arithmetic. IEEE Trans. Comput. **66**, 1281–1292 (2017). https://doi.org/10.1109/TC.2017.2690633

23. Fieker, C., Hart, W., Hofmann, T., Johansson, F.: Nemo/hecke. Proceedings of the 2017 ACM on International Symposium on Symbolic and Algebraic Computation - ISSAC '17 (2017). https://doi.org/10.1145/3087604.3087611

24. Johansson, F.: Computing hypergeometric functions rigorously. ACM Trans. Math. Softw. **45**(3), (2019). https://doi.org/10.1145/3328732

25. Johansson, F.: Faster arbitrary-precision dot product and matrix multiplication. In: 2019 IEEE 26th Symposium on Computer Arithmetic (ARITH), pp. 15–22 (2019). https://doi.org/10.1109/ARITH.2019.00012

26. Litzkow, M.J., Livny, M., Mutka, M.W.: Condor-a hunter of idle workstations. Proceedings of the 8th International Conference of Distributed Computing Systems, 104–111 (1988)

27. Bender, C.M., Orszag, S.A.: Advanced mathematical methods for scientists and engineers - asymptotic methods and perturbation theory. Springer-Verlag, ISBN 0-387-98931-5 (1999)

28. Demmel, J., Koev, P.: The accurate and efficient solution of a totally positive generalized Vandermonde linear system. SIAM J. Matrix Anal. Appl. **27**(1), 142–152 (2005)

29. The PARI group: PARI/GP Version 2.11.2. Univ. Bordeaux (2019). The PARI group. available from http://pari.math.u-bordeaux.fr/

30. Lenstra, A.K., Lenstra, H.W., Lovász, L.: Factoring polynomials with rational coefficients. Math. Ann. (1982). https://doi.org/10.1007/BF01457454

31. Schreiber, M.: Star polygons (2011). https://demonstrations.wolfram.com/StarPolygons/. Accessed 23 Sept 2022

32. Laura, P.A.: On the determination of the natural frequency of 8 star-shaped membrane. J. Royal Aeronaut. Soc. Pol. **68**, 274–275 (1964)

33. Wagner, H.: Fundamental frequency of a star-shaped membrane. Zeitschrift Angewandte Mathematik und Mechanik. **51**(4), 325–326 (1971). https://doi.org/10.1002/zamm.19710510417
34. Guidotti, P., Lambers, J.V.: Eigenvalue characterization and computation for the Laplacian on general 2-d domains. Numer. Funct. Anal. Optim. 29(5-6), 507-531 (2008). https://doi.org/10.1080/01630560802099233

35. Chen, C.S., Jiang, X., Chen, W., Yao, G.: Fast solution for solving the modified Helmholtz equation with the method of fundamental solutions. Communications in Computational Physics. 17(3), 867-886 (2015). https://doi.org/10.4208/cicp.181113.241014a

36. Kaiser, N.: Mean eigenvalues for simple, simply connected, compact lie groups. Journal of Physics A: Mathematical and General. 39(49), 15287-15298 (2006). https://doi.org/10.1088/0305-4470/39/49/013

37. Oikonomou, V.K.: Casimir energy for a regular polygon with Dirichlet boundaries (2010). arXiv:1012.5376

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.