On the Polygonal Faber-Krahn Inequality

It has been conjectured by P\'{o}lya and Szeg\"{o} seventy years ago that the planar set which minimizes the first eigenvalue of the Dirichlet-Laplace operator among polygons with $n$ sides and fixed area is the regular polygon. Despite its apparent simplicity, this result has only been proved for triangles and quadrilaterals. In this paper we prove that for each $n \ge 5$ the proof of the conjecture can be reduced to a finite number of certified numerical computations. Moreover, the local minimality of the regular polygon can be reduced to a single numerical computation. For $n=5, 6,7, 8$ we perform this computation and certify the numerical approximation by finite elements, up to machine errors.

In their book of 1951, Pólya and Szegö have conjectured a polygonal version of this inequality (see [42, page 158]). Precisely, denote by P n the family of simple polygons with n sides in R 2 and for every n ≥ 3 consider the problem (2) min This is the most technical part of the paper. Second, we give a bound for the maximal possible diameter of the optimal polygon as well as for the minimal length edge and inradius, when its area is fixed. As a consequence, it remains to prove that all polygons with free vertices in a (computable) compact set K ⊆ R 2n−4 are not optimal. This can be done by performing a finite number of certified computations of first eigenvalues, areas and perimeters. Indeed, if a polygon has vertices in the compact set K and is not optimal, then either due to uniform estimates of the modulus of continuity of the eigenvalue and measure or to monotonicty of both these quantities to inclusions, non-optimality is certified in an open neigbourhood. A finite number of such (open) neigbourdhoods will cover K.
Le us detail our strategy.
Step 1. (Formal computation of the shape Hessian matrix). We interpret the first eigenvalue as a function depending on the coordinates of the vertices of the n-gon (obtaining a function defined on a subset of R 2n ) and choose an appropriate, equivalent scale invariant formulation for problem (2). Once the validity of the first order optimality condition on the regular polygon is established, we compute the analytic expression of the shape Hessian. For that purpose, we rely on the computations done by A. Laurain in [36] for the energy functional (we recall the corresponding result in Remark 7.6) and perform similar computations for the eigenvalue, following the same method. Taking perturbations of polygons with n sides in the second shape derivative, we obtain the Hessian matrix (of size 2n × 2n) for the eigenvalue having the vertex coordinates as variables.
Step 2. (Numerical proof of the positivity of the shape Hessian matrix for the regular polygon, for a given n). The shape Hessian matrix of the scale invariant functional has four eigenvalues equal to 0, corresponding to the rigid motions and homotheties of the polygon. We use interval arithmetics and explicit error estimates for the finite element approximation to certify the positivity of the other eigenvalues of the shape Hessian matrix for the regular polygon with n sides. For n = 5, 6, 7, 8 and a suitable choice of an appropriate discretization, we certify, up to machine errors appearing in the meshing, the assembly and the resolution of the linear systems in the finite element method, that the remaining 2n − 4 of the eigenvalues of the Hessian are strictly positive. A fully certified (including machine errors aspects) positivity of the eigenvalues of the shape Hessian matrix is enough to prove the local minimality of the regular polygon, provided one knows that the coefficients of the matrix are continuous for small geometric perturbations of the regular polygon. This type of stability result is necessary to establish that the non zero eigenvalues remain positive in small neighborhood of the regular polygon. This is discussed in Step 3, below. By strict convexity, the regular polygon will be a minimizer in this neighborhood.
Step 3. (Quantitative stability of the shape Hessian matrix coefficients). Our objective is to identify a computable neighborhood of the regular polygon where the eigenvalues of a (2n − 4) × (2n − 4) submatrix of the shape Hessian matrix remains positive. The most technical part is to give analytic, computable, estimates of the variation of the coefficients of the Hessian matrix, for perturbations of the regular polygon. The difficulty comes from the fact that the expression of the coefficients involve the solutions of some (degenerate) elliptic PDEs with data in H −1+γ , depending on traces of the gradient of the eigenfunctions on segments. The analysis requires quantitative estimates of the perturbation of the eigenfunction in H 2 which relies, via Gagliardo-Nirenberg interpolation inequalities, on control of their norm in H 2+s . These estimates show that the unique, certified, computation of the Hessian matrix on the regular polygon is enough to obtain local minimality on a computable neignbourhood!
Step 4. (Analytic estimates of the maximal and minimal edge lengths of an optimal polygon). We give a computable estimate of the maximal diameter of the optimal polygon, provided its area is fixed. The estimate is inductively obtained for n ≥ 5: if the diameter of an n-gon exceeds some (computable) value, then its eigenvalue is close to the one associated to a polygon with n − 1 sides, so it can not be optimal in the class P n . Here we use surgery techniques inspired from [11], but face the difficulty of keeping constant the number of sides within the surgery procedure. As well, we give an analytic estimate for the minimal length of an edge and of the minimal inradius.
Step 5. (Formal proof of the conjecture). We show how to give an inductive formal proof of the conjecture reducing it to a finite number of (certified) numerical computations for each value of n. Up to this point we have computed, for the scale invariant functional, a neighborhood of the regular polygon where its minimality occurs and we have computed the maximal and minimal legnths of edges of an optimal polygon at prescribed area. Therefore, we are able to reduce the study of the conjecture to a family of polygons with vertices belonging to a compact set. Any certified evaluation of the eigenvalue/area of such a polygon showing non optimality, would readily produce a small neigbourhood of non optimal polygons, the size of the neigbourhood being uniform and analitically computed. Monotonicity with respect to inclusions of both the eigenvalue and the area may be very useful from a practical point of view, but not necessary for a theoretical argument. Finally we get a ball covering of a compact set which with known diameter, by balls of uniform size. This means that one can prove the conjecture after a finite number of numerical computations. We shall describe this procedure in Section 7.
This type of numerical procedure has successfully been used in [9] (to which we refer for a detailed description), for a different problem involving the same variational quantities but with only two degrees of freedom. The arguments transfer directly to our problem.
Although we prove that for a specific n the proof of the conjecture is reduced to a finite number computations, it is not our purpose to perform these computations, for two reasons. On the one hand, all constants that we prove to exist should be optimized and effectively computed. On the other hand, even for n = 5, in our procedure the number of degrees of freedom for the free vertices is 6 (see Section 7). This demands huge computational capacities. In other words, before any computational tentative, some further, deep, analysis should be performed to dramatically reduce the size of the computational tasks.
The structure of the paper is the following. Section 2 is devoted to the computation of the shape Hessian of the area and first eigenvalue functionals by a distributed formula. In particular, on polygons, we give the expression of the Hessian matrix of the eigenvalue as function of vertices coordinates. This section is inspired by the recent work of Laurain [36] for the energy functional. Section 3 contains a quantitative geometric stability result of the coefficients of the Hessian matrix with respect to vertex perturbations. This part is the key for the proof of the local minimality of the regular polygon and allows to estimate the size of the neighborhood of the regular polygon where minimality occurs. Sections 4 and 5 are devoted to the analysis of the shape Hessian matrix coefficients and to estimates regarding their numerical approximation. Section 6 contains certified computation of the eigenvalues of the shape Hessian matrix on the regular polygon, justifying, up to machine errors, its local minimality for n = 5, 6, 7, 8. In Section 7 we give an estimate of the maximal diameter of an optimal polygon and show how the proof of the conjecture reduces, for every n ∈ N, to a finite number of numerical computations. As well, we make short comment about the polygonal Saint-Venant inequality for the torsional rigidity, which can be analyzed in a similar way.

First and Second order shape derivatives
In this section we analyze the first and the second order shape derivatives of the first Dirichlet eigenvalue, for both general domains and for polygons. This section follows the strategy developed by Laurain in [36] for energy functionals (see Remark 7.6 for a brief summary of the corresponding results). Many proofs are very similar and we shall not reproduce them, referring to [36], whenever necessary. Nevertheless, the formulae for the eigenvalues are different, so that we shall detail them. The ultimate objective for polygons is to get an expression of the Hessian in distributed form involving sums over the two dimensional domains and remove any boundary integral expression. This is somehow contrary to what usually one does in shape optimization, the main motivation being that the distributed expression of the second shape derivative requires less regularity hypotheses than the boundary expressions. This is particularly useful for polygons. Finally, when restricted to the class of polygons with n sides, we shall describe the shape Hessian of the eigenvalue by a square symmetric matrix of size 2n × 2n.
In the literature one can find detailed descriptions of the shape gradients and shape Hessians of the eigenvalue on a smooth set (see for instance [26,27,34,13]). The case of polygons is more delicate, since the boundary expression of the shape Hessian fails to have sense, due to the lack of regularity of the boundary.
In order to simplify the reading and the interpretation of potential connections between the results of this section and [36], we use the same notations and, when the computations are similar, we prefer not to reproduce them and refer precisely to various sections in [36].
2.1. General domains. For vectors a, b ∈ R d and matrices S, T ∈ R d×d define the following: • Id denotes the identity matrix • a ⊗ b is the second order tensor of two vectors (a ⊗ b) ij = a i b j • a b = 1 2 (a ⊗ b + b ⊗ a) is the symmetric outer product. • a · b is the usual scalar product • S : T = n i,j=1 S ij T ij is the matrix dot product. It is immediate to notice that (a ⊗ b)c = (c · b)a and S : (a ⊗ b) = a · Sb.
As discussed in [36,Section 9.1], when computing second order shape derivatives, several approaches are possible. The one detailed in [36] uses the Eulerian derivative in order to compute the Fréchet derivative. However, the Eulerian derivative requires more regularity on one of the perturbation fields than W 1,∞ (R 2 , R 2 ), while perturbations of polygons are precisely in For a given vector field θ ∈ W 1,∞ (R d , R d ) consider the domain Ω θ = (I + θ)(Ω). It is well known that for θ W 1,∞ < 1 this transformation is an invertible diffeomorphism. In the following, when dealing with boundary value problems, we use subscripts to denote functions ϕ θ ∈ H 1 0 (Ω θ ) and superscripts to denote the functions ϕ θ = ϕ θ • (I + θ) ∈ H 1 0 (Ω). The objective in the following is to have distributed expressions which require less regularity than the generally well known boundary expressions for the shape derivative of the eigenvalue ( [27], [26]). Following the strategy of Laurain for the energy functional, we state below analogue results for the first and second Fréchet shape derivatives for the simple eigenvalues of the Dirichlet-Laplace problem (1). While some of these facts are standard (for instance the expression of the first derivative), the expression of the Fréchet second derivative and the matrix representation in the case of polygons seem to be new.
With these notations we are ready to state the following result.
Let λ be a simple eigenvalue of the Dirichlet Laplacian and u an associated L 2 -normalized eigenfunction. Then (i) The distributed shape derivative of λ is given by If, in addition, u ∈ H 2 (Ω), the corresponding boundary expression is (ii) The second order distributed Fréchet derivative is given by whereu(θ) andu(ξ) are the material derivatives in directions θ, ξ, respectively.
The first point is standard and may be found in many classical references, for instance [26]. Some formulae for the second derivative are also available in the literature, see [26], [27].The key point is that the distributed expression shown above is valid for Lipschitz domains and Lipschitz perturbations. Moreover, being written in symmetric form its expression helps in the computation of the Hessian matrix in the case of polygons.
Proof of Theorem 2.1. The first application of formula (6) is the expression of the first shape derivative. This computation is a classical result, but we present it here for the sake of completeness, since it illustrates well the techniques used when computing shape derivatives. Take v = u in (6) and note that, since u is the eigenfunction associated to λ(Ω), Using Ω u 2 dx = 1, we obtain A direct computation leads to Now we choose ξ ∈ W 1,∞ and we redo the same procedure to differentiate the first shape derivative (8). Denote Ω ξ = (I + ξ)(Ω) and suppose that ξ is small enough such that λ(Ω ξ ) is still a simple eigenvalue. Denote with u ξ ∈ H 1 0 (Ω ξ ) the eigenfunction associated to the simple eigenvalue λ(Ω ξ ). We have We also have the following elementary computation: As before, via a change of variables we write λ (Ω ξ )(θ) as an integral on Ω defining u ξ = u ξ • (I + ξ) ∈ H 1 0 (Ω). Using (4) and performing a change of variables, we obtain Now we are ready to compute the second Fréchet derivative of λ(Ω) by differentiating the previous expression w.r.t. ξ at 0 and denoting the derivative of u ξ at 0 byu(ξ). We use the product rule, differentiating the first term, the term D((I + ξ) −1 ) and finally det(I + Dξ). In particular, we have • det(I + Dζ) ζ (0)(ξ) = div(ξ). We obtain the following initial formula for the second shape derivative: Following [36, pag 25], we have −4(∇u(ξ) ∇u) : Dθ + 2∇u(ξ) · ∇u div θ = 2A (0)(θ)∇u · ∇u(ξ) and the material derivative (6) gives The derivative of the normalization condition gives We also have S λ 1 : DθDξ = (|∇u| 2 − λu 2 )Dθ T : Dξ − 2DθDξ∇u · ∇u, since tr(DθDξ) = Dθ T : Dξ. Combining all these expressions we obtain We have 2Dξ T ∇u ∇u : Dθ = (DξDθ + DξDθ T )∇u · ∇u.
Which gives This finishes the proof of the theorem.

2.2.
Polygons. In order to exploit the expression of Theorem 2.1 in the case when Ω is a polygon, we follow again the strategy of Laurain [36] to extend a geometric perturbation of vertices to a global perturbation of the polygon.

Vertex perturbation versus global perturbation.
Suppose Ω is a n-gon. Starting from a perturbation of the vertices, the perturbation field θ ∈ W 1,∞ (R 2 ) will be built as follows.
Denote the vertices of the polygon by a i ∈ R 2 , i = 0, ..., n − 1 and for each vertex consider the vector perturbation θ i ∈ R 2 , i = 0, ..., n − 1. Whenever necessary, we suppose that the indices are considered modulo n. Consider a triangulation T of Ω such that the edges of the polygon are complete edges of some triangles in this triangulation. Moreover, consider the following globally Lipschitz functions ϕ i for 0 ≤ i ≤ n − 1 that are piecewise affine on each triangle of T and satisfy Several choices are possible, as the two examples of Figure 1 show, their extension outside the polygon being irrelevant. Then, we build a global perturbation of R 2 given by Gradient and Hessian of the area functional. The shape derivatives for the area functional are classical and are widely studied in the literature (see [26], [36], etc.). The expression of the shape derivative of the area is However, in the particular case of n-gons the situation is much simpler, since explicit formulae exist in terms of the coordinates of the vertices of the polygon. For a non degenerate polygon whose coordinates of the vertices are denoted by (x i , y i ) and whose edges are oriented in the counter-clockwise order the area is given by The coordinates are regrouped in the vector by concatenating the coordinates of the vertices a i (13) x = (a 0 , a 1 , ..., a n−1 ) = (x 0 , y 0 , ..., x n−1 , y n−1 ) ∈ R 2n , which will always be the case in the following, when parametrizing polygons. The gradient of the area in terms of the coordinates verifies: We denote by R c,α the rotation around c ∈ R 2 with angle α (in the trigonometric sense), hence the gradient of the area has the geometric expression This is natural, since the area of the polygon when moving a vertex a i only varies when moving the vertex a i in the normal direction to the closest diagonal. Another expression of the gradient of the area, using the functions ϕ i defined earlier, can be found following the results of [36] and is given by (15) ∇A Since the expression of the gradient of the area is linear in terms of the coordinates, the Hessian matrix of the area of the polygon is the constant 2n × 2n block matrix where the non-zero 2 × 2 blocks are given by Figure 2. Boundary perturbation induced when perturbing a vertex Following the results in [36] we find that the formula for the Hessian of the area in terms of the functions ϕ i can also be expressed using the following block structure In particular the Hessian of the area can be written as a tensorial product (Kronecker product) between the matrices  Therefore, the corresponding eigenvalues and eigenvectors can be found explicitly.
Gradient of the eigenvalue. Below we compute the gradient of the eigenvalue (1) as function of the vertices, i.e. the partial derivatives of these functionals with respect to the coordinates of the vertices of the polygons. The expression of these gradients can be used to prove that the regular polygon is a critical point under an area constraint and are useful for numerical computations.
The expression of the gradient of the eigenvalue with respect to the coordinates is a consequence of the shape derivative formulae recalled in the previous section. It is enough to use the distributed expression of the shape derivative, valid in general, with the perturbation field θ introduced in (11). An example is given in Figure 2 for θ = θ i ϕ i . The proof is similar to the case of the torsion energy [36]. We choose to detail here only the boundary expression, with a slightly different argument than the one used in [36].
Theorem 2.2. The gradient of a simple Dirichlet-Laplace eigenvalue (1) when Ω is a polygon with coordinates x as in (13) is given by where n is the outer unit normal vector.
Notice that the boundary expression is always valid, even though the eigenfunction itself does not belong to H 2 (Ω). This is a consequence of the fact that in an arbitrarry polygon (typically non convex), the eigenfunction enjoys a local H 2+δ regularity far from the corners, while at corners the singular part has a very specific structure, albeit good enough to make the boundary expression of the gradient valid. We recall from [7] that where u reg ∈ H 2+δ (Ω) for some δ > 0 and where C i are constants, ω i are the angles, ψ i is cutoff function equal to 1 in a neighborhood of the vertex a i and (r, θ) are the polar coordinates around the angle i.
Proof. The expression ∇λ(x) = Ω S λ 1 ∇ϕ i dx i=0,...,n−1 is valid. It remains to prove the equal- First, note that the gradient of u is point-wise defined on ∂Ω, except at the vertices, in a classical way. We fix a vertex i and define Since u| Ωε ∈ H 2 (Ω ε ), a direct computation shows that div S λ 1 = 0 on Ω ε , the divergence being applied on lines. Moreover, since u = 0 on ∂Ω, the gradient ∇u is colinear with the normal vector n on ∂Ω. In particular, (∇u ⊗ ∇u)n = (n · ∇u)∇u = |∇u| 2 n. As a consequence S λ 1 n = −|∇u| 2 n. Therefore we obtain We conclude by noticing that which is a consequence of the decomposition u = u reg + u sing . We know that u reg ∈ H 2+δ (Ω) and H 2+δ (Ω) is embedded in W 1,∞ (Ω), so that the gradient of u reg is bounded. At the same time, |∇u sing | ≤ Cr π ω i −1 for some constant C independent on ε. Both these observations lead to Γε |∇u reg | 2 + |∇u sing | 2 → 0, for ε → 0.
To conclude notice that Remark 2.3. It is possible to note that the integrals which come into play in the boundary expression of the gradient only need to be computed on two adjacent sides to vertex a i , which gives In the following we make the convention that the Jacobian matrix of a vector function contains gradients of the components on every line. Hessian matrix of the eigenvalue. Following the notation of [36], we introduce the functions U i ∈ H 1 0 (Ω, R 2 ), i = 0, ..., n − 1 such thatu(θ) = n−1 i=0 θ i · U i . Using (6) we get the set of two PDEs: , for every v ∈ H 1 0 (Ω). The normalization condition (7) gives (19) Ω (2uU i + u 2 ∇ϕ i ) dx = 0, so that the system of equations (18) - (19) has a unique solution U i .
Proof of Theorem 2.4: The proof of this result, is computational in nature and is inspired by [36,Proposition 14]. To obtain the Hessian matrix we use the formula for K λ given in Theorem 2.1 for the Fréchet second shape derivative. There are several terms, already computed in [36, Appendix A], which also appear in the formula for the eigenvalue. We only present in detail the terms which are different. We point out that in order to obtain directly the Hessian matrix, the 2 × 2 blocks should be multiplied by the variables ξ j below, which gives transposed 2 × 2 blocks compared to [36].
The first term is straightforward The second term is treated in [36] (term L 3 , pag. 38): The third term is similar to the term L 4 treated in [36] (pag. 39): The fourth term treated in [36] (L 5 pag. 39): 2 The fifth term is new and will be computed below. Note that under the conventions θ = n−1 i=0 θ i ϕ i and ξ = n−1 i=0 ξ i ϕ i (see (10)- (11) for the definition of ϕ i ) we have: Regrouping all the above results finishes the proof of the theorem.
Remark 2.5. It is worth to notice that the matrix N λ obtained in Theorem 2.4 and the corresponding matrix obtained by Laurain in [36,Proposition 14] have similar structures (see Remark 7.6). Moreover, the results resemble the structure of the tensor S λ 1 corresponding to the first shape derivative in distributed form. The matrix N λ has an additional term coming from the fact that the eigenvalue λ(Ω) is already present in S λ 1 , and its derivative appears when computing the second shape derivative.
Remark 2.6. It can be noted that the Hessian matrix found in (20) does not depend on the normalization condition (19). It is more convenient in the following to suppose that the functions U i are normalized with the following condition where u is the eigenfunction associated to the simple eigenvalue λ(Ω) of the Dirichlet-Laplacian.
General properties of the Hessian matrix. The formulas for the gradient and the Hessian matrix obtained previously do not depend on the choice of the perturbation given in (11). As illustrated in Figure 1 multiple choices for the triangulations defining the functions ϕ i are possible. In particular: • when the triangulation contains no inner vertices then n i=1 ϕ i = 1, which implies that n i=1 ∇ϕ i = 0. • for the regular polygon, considering a triangulation with an additional vertex at the center of the polygon provides additional symmetry properties. In the following we will switch between the two choices above in order to obtain further properties of the gradient and the Hessian matrix. In the following, define the two vectors t x = (1, 0, 1, 0, ..., 1, 0) and t y = (0, 1, 0, 1, ..., 0, 1) ∈ R 2n . Proposition 2.7. 1. The sum of the components on of the gradient ∇λ(x) on odd and even positions, respectively is zero. Equivalently we have ∇λ(x) · t x = ∇λ(x) · t y = 0.
2. The vectors t x , t y are eigenvectors of the matrix N λ defined in (20).
Proof. Let us note that by choosing ϕ i on a triangulation with no interior vertices we have n i=1 ∇ϕ i = 0. This already gives an answer to the first point above since For the second point, let us note that with the same choice of the functions ϕ i the solutions U i of (18) with the normalization condition (21) verify n i=1 U i = 0 since the sum of the right hand sides in (18) is equal to zero. It is now straightforward to see that N λ t x = N λ t y = 0 which implies that the vectors t x , t y are eigenvectors of N λ corresponding to the zero eigenvalue.
Formula (20) respects the structure of the second shape derivative. It is possible to simplify the formula using the definition of S λ 1 and the property (a ⊗ b)(c ⊗ d) = (b · c)(a ⊗ d). Regrouping terms we obtain It is immediate to see that Therefore, the expression of the Hessian matrix simplifies to From this point on, in the rest of the paper, we concentrate on the case of the first eigenvalue of the regular polygon and we further simplify the expression of the Hessian. By uniqueness arguments the first eigenfunction u of the Dirichlet Laplace operator on the regular polygon has the same symmetries as the regular polygon.
In the following suppose that ϕ i , 0 ≤ i ≤ n − 1 are associated to the particular triangulation T = (T k ) n−1 k=0 of the regular polygon made of congruent triangles with one vertex at the center (see Figure 1). Thus, the triangulation T also respects the symmetry of the regular polygon. The symmetry of the first eigenfunction implies that T k (|∇u 1 | 2 − λ 1 (Ω)u 2 1 ) dx = 0. Using this relation the gradient of λ 1 (Ω) on the regular polygon becomes Using the fact that ∇ϕ i ⊗ ∇ϕ j − ∇ϕ j ⊗ ∇ϕ i is piece-wise constant on every triangle T k , k = 0, ..., n − 1, we find that where B ij are the blocks of the Hessian of the area given in (17).
Recall that ∇ϕ i is piecewise constant on the triangles T k and by symmetry T k u 2 1 dx = 1/n for k = 0, ..., n − 1. Therefore T k u 2 is the area of the polygon having vertices at coordinates given by x, as recalled earlier. Therefore, the last term in N λ ij has the form Consider now the Hessian of the product λ 1 (x)A(x) and note that we have In this formula the last term of the Hessian of λ 1 (x) simplifies the tensorial products between the gradient of the area and the gradient of the eigenvalue. Following the previous computations we arrive at the following significant simplification for the Hessian of the product of the area and the eigenvalue.
Proposition 2.8. In the case where Ω is a regular n-gon and the triangulation T defining ϕ i is symmetric the Hessian matrix of λ 1 (Ω)|Ω| = A(x)λ 1 (x) in terms of the coordinates of the polygon has the 2 × 2 blocks M λ ij , 0 ≤ i, j ≤ n − 1 given by The simplified formula (23) for the Hessian of the product of the area and the first eigenvalue has three terms: • The first one is related to the decomposition U i of the material derivatives given in (18). Furthermore, the terms are related to the bilinear form from the variational formulations of U i , which will be essential in improving the estimates in the numerical simulations. This part of the Hessian is negative definite. • The second term is related to the Hessian of the area given in (17). The associated blocks are non-zero only when |i − j| = 1 (modulo n). This part has both positive and negative eigenvalues. • The third term involves only the first eigenfunction u 1 and the functions ϕ i defined in (10). The associated blocks are non-zero only when |i − j| ≤ 1. This part of the Hessian is positive definite. Although the expression of the Hessian given in (23) is explicit, its positive definiteness is not obvious. The analysis of the eigenvalues of this matrix is continued in Section 4.

Geometric stability of the shape Hessian matrix
In this section we shall perform both a qualitative and quantitative analysis of the behavior of the coefficients of the Hessian matrix for local perturbations of the vertices of the regular polygon P n inscribed in the unit circle with one vertex at (1, 0). Some of the results would extend naturally either to perturbations of general convex polygons or even to more general sets. Nevertheless, we focus on the perturbation of the regular n-gon and we shall not search generality. The two main technical aspects of this section are described below.
• Continuity of the Hessian matrix coefficients for the geometric perturbation.
We prove the continuity of the shape Hessian matrix for a perturbation of the regular polygon. This question is itself non trivial because of the weak regularity of the right hand sides in the equations satisfied by the solutions U i of (18). Stability results for the eigenfunctions in H 2 are required, whereas the classically known stability based on γ-convergence holds in H 1 . The continuity of the coefficients will readily give the local minimality of the regular polygon provided the positive definiteness of the Hessian matrix is known on the regular polygon only. • Estimate of the modulus of continuity of the coefficients for the geometric perturbation. This information is crucial to formally reduce the proof of the conjecture to a finite number of numerical computations. We compute the modulus of continuity of the coefficients, i.e. we find estimates of the variation of all coefficients of the Hessian matrix in terms of some power of Hausdorff distance between the perturbed polygon and the regular polygon. In other words, for every δ > 0 we identify a value ε > 0 such that all the coefficients of the Hessian matrix computed on polygons with n sides in an ε-neigbourhood of P n stay in a δ-neighborhood of the coefficients of the Hessian matrix associated to P n . We split this section in three subsections, going from basic estimates for the variations of the eigenvalues and eigenfunctions to the estimates of the variation of the matrix coefficients. This last point is more delicate as it involves solutions of (18)- (21) with variable, singular, right hand sides that are not in L 2 .
Throughout this section, we denote by C, θ two positive constants which may change from line to line. The tracking of those constants is possible but, since we will not perform here numerical computations of an effective neighborhood of minimality, this is not immediately useful. Consequently, in order to avoid heavy calculations we choose to prove only the existence of those constants. In particular, we are not aimed here to optimize the constants, which in case of certified numerical computations of the neighborhood would be a priority.
In the particular case in which f = 1, we denote w Ω the the solution of (24), and call it torsion function. The torsion function is the unique minimizer of the torsion energy, Let now Ω α , α ∈ {a, b} be two such domains and denote by v α the solution of (24) on Ω α for the right hand side f α and by u 1,α the L 2 -normalized, non-negative eigenfunctions on Ω α corresponding to the first eigenvalues λ 1,α , respectively. We denote by d H the Hausdorff distance.
In a first step, we seek estimates of the form Above, all functions u 1,α , v α are assumed to be extended by 0 on the complement of their definition domain, this extension being suitable for H 1 -estimates. By abuse of notation, the extensions by 0 are still denoted with the same symbols. The literature is quite rich for such type of H 1 -estimates, like (25) and (27). For instance, Savaré and Schimperna [45] give estimates for solutions of (24) in the class of sets satisfying a uniform cone condition while Burenkov and Lamberti [5], Feleqi [18] discuss the eigenfunctions. Concerning (26), we refer to [41] (see as well Section 7) for sharp estimates with power ϑ = 1 2 and controlled constant. Let us point out a relevant fact, which becomes important as soon as we search to identify all the constants in (25)- (27). The results referred above occur in the class of domains satisfying a uniform cone condition, while our setting is much more regular: we locally perturb the regular n-gon, always obtaining a convex n-gon. This regular behavior will be exploited in the next subsection to get estimates in higher order norm even in the case of singular right hand sides and it dramatically simplifies the proofs of the H 1 -estimates.
Below we shall only recall some results without proofs. The interested reader could easily recover the estimates in our regular setting in a more direct way. Assume that Ω a , Ω b ⊆ R 2 satisfy a uniform (ρ, ε)-cone condition (see [45,Definition 2.6]). [45]). If f a = f b := f , there exists a constant depending only on the diameters such that Note that the first two inequalites require f ∈ L 2 (R 2 ). The result recalled in Proposition 3.1 together with the Poincaré inequality readily gives inequality (25). Note as well that the Poincaré constants on the two domains equal the first Dirichlet eigenvalues.
For a small perturbation of the regular n-gon, the values of ρ and ϑ can be computed explicitly. However, in this last case a more direct proof of the inequalities can be obtained as a consequence of the uniform bound of the H 2 norms of the solutions with an explicit value (maybe not optimal) of the constant C.
Concerning the estimates (26) and (27), we refer to the papers of Feleqi [18] and Burenkov and Lamberti [12]. Those esimates being less explicit, we give below a slef contained argument which takes advantage of the convexity of the sets.
For now, assume that Ω a and Ω b are convex, in which case the level sets of the torsion function and of the first eigenfunctions are convex. Moreover, v α and the eigenfunction u 1,α belong to H 2 (Ω α ) as we shall recall in the next subsection. We recall a first regularity result in the class of convex sets, due to Grisvard [ For a n-gon which is a small perturbation of the regular n-gon P n , this inequality gives uniform bounds for the H 2 -norms of the normalized eigenfunctions and of some H 2 extensions in R 2 . The bounds in L ∞ are standard and the convexity of the polygon together with the barrier method provides L ∞ estimates for the gradients.
Denotingṽ the solution of (24) inΩ, we have We notice that the function v α −ṽ is harmonic onΩ, so its maximum onΩ is attained on ∂Ω, whereṽ vanishes. Since ∂Ω lies in a neighborhood of ∂Ω α , denoting However, for every In order to bound ∇w α ∞ we take advantage that the level sets of w α are convex and so we have a barrier given by the width. Indeed, in every point x of the the level set, we can find an infinite strip containing the level set and having one boundary line passing through x. Using the classical barrier method gives |∇w α (x)| ≤ W x /2, where W x is the width. This implies Adding the estimates for v a and v b leads to the conclusion.
Perturbations of the regular polygon. For n ≥ 5 we denote P n = a * 0 a * 1 . . . a * n−1 the regular polygon with n sides inscribed in the unit circle with a * 0 = (1, 0). We denote R n , r n the radii of the circumscribed, inscribed circles for P n and l n the length of an edge. Denote the area of P n by A n . The angles are equal to n−2 n π. An easy computation leads to R n = 1, r n = cos π n , l n = 2 sin π n , A n = n sin(2π/n)/2.
Let P denote generically a perturbation of P n , i.e. polygon a 0 a 1 . . . a n−1 with n sides such that the angles of the perturbed polygon do not exceed We can represent both the boundaries of P n and P using the same n charts given by the graphs of the boundaries ∂P n , ∂P over the segments In each chart, the function representing the boundary of the polygons is piecewise affine with two slopes not exceeding tan( π n + arctan( 1 4 sin π n )). For n ≥ 5 an upper bound for this quantity is 0.73.
We denote λ k , λ * k the k-th eigenvalues and u k and u * k the corresponding normalized eigenfunctions on P , P n , respectively.

Proposition 3.4. Under the previous hypotheses
Proof. The inclusions rn rn+ε P ⊆ P n ⊆ rn rn−ε P , imply that r N +ε We introduce the problem

Using Lemma 3.3 and the Poincaré inequality we have
.
Using the orthonormal Hilbert basis of eigenfunctions in H 1 0 (P ) we consider the decomposition We have On the other hand, which, after elementary computations leads to Finally, By summation, the inequality follows.
Remark 3.5. In order to complete the estimates we recall that in simply connected domains Grebenkov [21,Formula (6.22)]). We also recall from [3] that λ 2 λ 1 ≤ j 2 1,1 /j 2 0,1 , where j 0,1 , j 1,1 denote the first positive zero of the Bessel functions J 0 , J 1 and that λ 2 − λ 1 ≥ 3π 2 diam 2 (P ) from [1]. As well, by inclusion and homogeneity, λ 1 (P n ⊕ B ε ) ≥ 1 1+ε 2 λ * 1 . We can also give a direct estimate for ψ − u 1 2 . Indeed, Proposition 3.6. There exists a constant C > 0 such that for all 0 < ε < ε 0 The first inequality is a consequence of the barrier method. The diameter and the inner ball control the size of the eigenvalue and of the L ∞ norm of the the eigenfunctions, themself being controlled by ε 0 .
The second inequality is a consequence the Gagliardo-Nirenberg inequality (see for instance [43])

3.2.
Uniform H 2+s regularity of the eigenfunctions. In this section we recall some finer estimates of the regularity of the solutions v α of (24) in polygons which are small perturbations of the regular polygon. However, we need more regularity than H 2 in order to quantify the variation of the shape Hessian coefficients. These finer regularity results take full advantage from the very specific convex, polygonal geometry of the domains, size of angles and number of local charts of the boundary. We refer the reader to [15] for detailed analysis of the regularity in polygonal domains.
Lemma 3.7. Let P be a perturbation of the regular polygon P n as above. Let 0 < γ ≤ π ω 0 . Then, for every f ∈ H −1+γ (P ) the solution of (24) The constant C depends on γ but it is independent on f and P .
Above, the independence on P comes precisely from the very specific perturbation we consider, which keeps constant the charts and controls the angles. Let us denote s 0 = π ω 0 − 1 > 0 and let 0 ≤ s ≤ s 0 .

Corollary 3.8. Under the previous hypotheses and notations we have
with C depends on s but is independent on the perturbation.
Proof. This is a consequence of Lemma 3.7 and of the fact that the right hand sides λ 1 u 1 of the equations solved by the eigenfunctions have an H 1 -norm equal to λ 1 (1 + λ 1 ) which is uniformly bounded in the class of perturbations we consider.
One has to pay particular attention to the extension of u 1 on the complement of P . As far as we are concerned with L p , H 1 properties of the extension, performing an extension by 0 on R 2 \ P is enough. Neverhtless, such an extension does not belong to H 2 , H 2+s , so we can not compare the extensions of u 1 and u * 1 in those norms. Two choices can be done in order to compare solutions on different polygons in H 2 . Either we extend them in H 2 and compare their extensions, or we locally compare on compact sets included in both domains. Below, we choose to compare their extensions. The extensions we seek rely on the Stein universal extension operator (see [47] and [33,29]). We recall the following from from [47]. Proposition 3.9. Assuming P is a perturbation of the regular polygon as above, there exists an extension operator where the constant C above depends on q but not on P .
Remark 3.10. We point out that the extension of Stein relies mainly on the construction of a smoothed distance function. The choice of this function is not unique. Stein proposed a construction based on partition of the complement of P on squares belonging to the union of latices (2 −k Z 2 ) k∈Z . In the sequel we shall use this argument and the freedom to build the smoothed distance function in order to be able to compare the extension operators on P and P n . Using a cut off function, we will assume that all extensions E P (u) vanish outside the ball B 2 .
Proposition 3.11. There exists C > 0, ϑ ∈ (0, 1) such that for every u ∈ H 2+s (R 2 ) The key use of this result is related to the possible extensions of an eigenfunction outisde P . Indeed, from Proposition 3.4 we control the norm u 1 − u * 1 H 1 (R 2 ) . However, this is true for the extensions by 0 of the eigenfunctions not for the extensions given by the Stein operator. Proposition 3.11 together with Proposition 3.9 imply that we can control the norm of the difference in H 2 for the Stein extensions provided we control the norm in L 2 . This is a consequence of the following Lemma.
Lemma 3.12. By E Pn we denote a (suitably chosen) Stein extension operator associated to P n . There exists a constant C such that for every perturbation P as above there exists a Stein extension operator E P satisfying . Proof. We rely on the construction of the operator by Stein using the averaging method (see [47,Theorem 5,page 181]). The difficulty is that we deal with extension operators corresponding to different domains and applied to different functions. We want to prove that the extended functions are close in L ∞ provided that the non extended functions are close in L ∞ . Since each one is extended with its own operator, we have to detail the construction of the operators in order to be able to perform the comparison.
Step 1. Localization. Since the boundary of P is described in the same charts as the boundary of the regular polygon, we use the explicit formula of the extension operator. We refer the reader to [47, Theorem 5, page 181] (see also [33,29]), where the explicit construction is given.
There exists a smooth partition of unity consisting on n + 2 functions (ψ j ) j=0,...,n+1 such that for every vertex a j of P n there exists one function ψ j supported in B(a j , 3 4 l n ), one of the functions is supported in Int(P n ) and one is supported in Int(R 2 \ P n ). In view of the smallness of the perturbation P of the regular polygon, we can keep the same n charts to describe the boundary of ∂P and use the same partition of unity as above, for the regular polygon. The maps of the charts are built in a uniform way as piecewise affine functions having two controlled slopes.
Moreover, instead of extending u 1 , u * 1 we shall extend each function u 1 ψ j , u * 1 ψ j relying on the special construction given by Stein in [47, Theorem 5, page 181], which takes advantage from the specific graph structure of the boundary. Finally, we use the generic comparison Step 2. Construction of the smoothed distance functions. The expression of the Stein extension operator is explicit and relies on regularization of the distance functions to P, P n respectively, say ∆ P , ∆ Pn . The construction of these functions is quite delicate and we refer the reader to [47, Theorem 2, page 171] for all the details. We have ∆ P ∈ C ∞ (R 2 \ P ), satisfying In its construction, Stein gives a precise formula for ∆ P , namely where Q k consists in a suitable partition of R 2 \ P in squares and φ k are C ∞ functions equal to 1 on Q k and vanishing outside a 9 8 -dilation of Q k by the center of Q k . The partition (Q k ) k is not arbitrary, the size of the squares being controlled by the distance of the square to the boundary of P .
Assume now that P is a perturbation of P n as above such that d H (∂P, ∂P n ) = ε. Then, Our aim is to slightly modify the construction of the partition (Q k ) for P such that at distance larger than 16ε from the boundary of P , the partition coincides with the one associated to P n . This will entail that if d(x, P ) > 128ε then ∆ P (x) = ∆ Pn (x). This is done as follows.
• We first set the family grids (2 −k Z 2 ) k∈Z in R 2 and choose a suitable partition for R 2 \ P n .
• We select out from this partition all the squares which intersect the set We use the Stein's method to fill the rest of the partition associated to P , namely to cover the open subset of R 2 \ P not yet covered by the selected partition. Finally, the construction of the functions φ k follows the same procedure as Stein. The only difference from the original Stein construction is only the alteration of the partition at distance larger than 16ε. In view of (36), properties (34)-(35) of ∆ P are preserved.
The main consequence of this construction is that if d(x, P ) > 128ε then ∆ P (x) = ∆ Pn (x).
Step 3. Comparison of the extensions. We recall that u 1 and u * 1 are uniformly Lipschitz in R 2 , as a consequence of Proposition 3.6. This plays a crucial role in estimate (33). Let us now recall from [47] how the Stein extension works. We shall simultaneously write the extension of u 1 with E P and the extension of u * 1 with E Pn . Suppose P, P n are above the graphs representing their boundaries on a segment [m j , M j ], which we suppose, without loss of generality, is contained in the horizontal coordinate axis.
Let c > 0 be a constant such that The extension operators are defined for x ∈ [m j , M j ] and y < φ j (x) and y < φ * j (x), respectively, by To complete the estimate, we evaluate both E P (u 1 )(x, y) and E Pn (u * 1 )(x, y)| for (x, y) lying at distance not larger than 130ε from the boundary of P n . Here we take advantage from the fact that there exists C, independent on P (see [47,Theorem 5, page 181]) such that . Since u 1 , u * 1 vanish on ∂P, ∂P n , respectively, we get that for (x, y) as above we have . This last inequality concludes the proof.
As a consequence of the Proposition 3.11 and Lemma 3.12, together with the uniform boundedness of the support of the extended functions, we get the following. Corollary 3.13. There exist constants C and ϑ ∈ (0, 1) independent on the perturbation, such that Estimates of the Hessian coefficients along the perturbation. In the sequel we collect some L ∞ -estimates, necessary for estimates of the coefficients of the Hessian matrix. Let ϕ * : T * → R, ϕ : T → R be the functions defined in (10) (the second kind, in Figure 1). We The last inequality takes advantage from the previous one and from the fact that the level sets are convex, via the barrier method.
Then, for every s ∈ (0, 1 2 ] there exists a constant C s depending only on s, such that Φ . In the last inequality, we used the classical trace inequality in H 1 (R 2 ) and the fractional trace inequality in H 1 2 +s (R 2 ) (see [48,Lemma 16.1]) together with the continuous embedding of H s (−1, 1) ⊆ L 2 (−1, 1).
Proof. We shall make an explicit computation. LetS 2 = [B 1 B 2 ] be the segment on the same line as S 2 such that its vertical projection on the horizontal axis is precisely S 1 . The A 1 B 1 ≤ ε and A 2 B 2 ≤ ε. We have the following estimates.
Let us introduce the projector Π 1 :S 2 (x, y) → (x, 0) ∈ S 1 . Then, Adding all the previous estimates, we conclude the lemma.
We turn our attention to U i , the solution of (18) - (21) in P . Recall that the expression of the coefficients of N ij in (20) does not change when a multiple of the eigenfunction u 1 is added to U i . In the following, whenever working with vectorial quantities, estimates are understood component by component.
We drop the index i and we formally write (18) and involves the following type of terms (possibly multiplied by geometric quantities) where S is an edge of T and n is the normal. Note that u 1 ∈ H 2+s (P ) and all these quantities are controlled for our perturbation, in a norm which is at least H −1+s .
Proof. One readily gets , which gives, using the Poincaré inequality in the orthogonal of u 1 , Taking into account the Andrews-Clutterbuck result [1] and the structure of f , Lemma 3.7 gives the conclusion.
In order to estimate we rely on the stability estimates for simultaneous domain and right hand side perturbations. Moreover, in view of the definitions of U * , U, we have to work in the orthogonal on u, u * , and use a correction term built by projection.
We have the following.
Lemma 3.17. There exist positive constants C, ϑ > 0, such that for every admissible perturbation Without restricting the generality we can assume that P ⊆ P n . Indeed, if this is not the case, we compare both U and U * with the solution U on the regular polygon (1 + ε)P n , which contains both P and P n .
We introduce the following auxiliary problem In the same time, both V and U * belong to H 1+s with controlled norm, so in particular they belong to W s 2 ,∞ with controlled norm. Using again the Gagliardo-Nirenberg inequality for the Stein extension of V, we get At the same time, −∆Ṽ − λṼ = λ * 1 U * + f * − λV := f in D (P ) and by straightforward computation Since we can use the stability result (39) to We can now conclude with the following.
Theorem 3.18. There exists C, ϑ > 0 such that for every polygon P ∈ P n satsifying ∀i = 1, . . . , n, Proof. The first inequality is a direct consquence of Lemma 3.17. The second one is a further consequence of the Weyl inequality on the stability of eigenvalues for perturbations of a symmetric matrix and on the equivalence of all norms over a finite dimensional space.
Remark 3. 19. The Hessian matrix of the area of the polygon is constant. As a direct consequence, a similar estimate holds for the Hessian matrix M λ of the scale invariant functional P → |P |λ 1 (P ).

Eigenvalues of the Hessian matrix for the regular polygon
We denote again P n = [a 0 a 2 ...a n−1 ] the regular polygon with n-sides, centered at the origin, with the vertex a 0 at the point (1, 0). As well, λ 1 := λ 1 (P n ) denotes its first eigenvalue and u 1 := u 1 (P n ) a positive, L 2 -normalized eigenfunction. We also use the notation θ = 2π/n.
As a consequence of the homogeneity of the eigenvalue to rescalings λ 1 (tP ) = 1 t 2 λ 1 (P ), the proposition below establishes the equivalence between the original problem (2) and some unconstrained versions. Its proof is standard and will not be recalled. For the convenience of the reader, we also collect below some well known facts.
Proposition 4.2. Let n ≥ 3. Then (1) The first eigenfunction on P n has the symmetry of the n-gon.
(2) • P n is a critical point for problem (L 1 ) above; • any regular n-gon is a critical point for problem (L 2 ) above (see Theorem 4.14); • the regular n-gon λ 1 (Pn) |Pn|c 1 4 P n is critical for problem (L 3 ) above. (3) If moreover any of the regular n-gons above is a local minima for its own problem, then all the others are local minima for their own problems.   [19] show that the equilateral one is the only possible critical point for the two functionals (first eigenvalue and torsional rigidity) studied here.
Proposition 4.5. Let P n be the regular polygon defined above. If the Hessian matrix M λ of P → |P |λ 1 (P ) evaluated at P n , given in (23), has 2n − 4 eigenvalues that are strictly positive then P n is a local minimum.
Proof. In the previous section in Theorem 3.18 it is shown that the coefficients of Hessian matrix are continuous for a local perturbation of the free vertices. Therefore, it would be enough to prove that the Hessian matrix associated to the free variables is positive definite. Fix the two consecutive vertices a n−2 , a n−1 and consider the associated matrix M which is the (2n − 4) × (2n − 4) principal submatrix of M λ obtained by removing the last four lines and columns. Then M is the Hessian matrix of the same functional, with the last four variables removed.
First of all, we observe that M λ has 4 zero eigenvalues which correspond to translations, scalings and rotations which leave the objective function invariant. In Propositions 4.6, 4.12 direct proofs are given showing that (41) t are indeed eigenvectors of M λ associated to the zero eigenvalue.
Suppose that M λ has 2n − 4 strictly positive eigenvalues (in addition to the four zero eigenvalues described above). The result stated in [30,Theorem 4.3.28] shows that the eigenvalues of M have lower bounds given by those of M λ , therefore they are non-negative. Suppose that M has a zero eigenvalues with an eigenvector ξ ∈ R 2n−4 . Completing ξ with zeros would give an eigenvector of M λ associated to the zero eigenvalue. This is impossible since taking the last four components of the eigenvectors in (41) gives four independent vectors in R 4 . Therefore M is positive definite implying that P n is indeed a local minimum for the functional P → λ 1 (P )|P |.
The remaining part of this section is dedicated to the computation of the eigenvalues of M λ . In particular, we show that the eigenvalues of M λ can be computed in terms of the first eigenfunction u 1 and the solutions (U 1 0 , U 2 0 ) of (18) with the normalization condition Pn U i 0 u 1 = 0, i = 1, 2. A numerical approach for proving that the matrix M λ has 2n − 4 eigenvalues that are strictly positive is provided in the next section. Proof. The proof is immediate, following the expression of M λ given in (23). Proposition 2.7 shows that t x and t y are in the kernel of the Hessian of the eigenvalue and are orthogonal to both the gradients of the eigenvalue and of the area. Moreover, they are also in the kernel of the area Hessian (16). Combining all these aspects finishes the proof.
The following result recalls the symmetry properties of u and U 1 0 , U 2 0 . For simplicity, we use the notation a(u, v) = Pn ∇u · ∇v − λ Pn uv. Proposition 4.7. The following holds.

The quantities
j → a(U 1 0 , U 1 j ), j → a(U 2 0 , U 2 j ) are even with respect to j (modulo n) and the quantities The proof is straightforward from the definitions. Change of basis. In order to deduce more information about the structure of the Hessian matrix it is useful to perform a change of basis so that for each vertex the basis directions correspond to the radial and tangential directions (see Figure 3). The Hessian matrix in the new basis is given by the formula H λ = P T M λ P where P = (P ij ) 1≤i,j≤n is a 2 × 2 block matrix with P jj = cos(j − 1)θ − sin(j − 1)θ sin(j − 1)θ cos(j − 1)θ . Of course, M λ and H λ have the same eigenvalues.
Moreover, the Hessian matrix H λ in this particular basis has an additional property. Indeed, it can be seen that in this basis the matrix does not change when a circular perturbation is  Figure 3. Change of basis to radial and tangential components (left). An example of symmetric triangulation defining ϕ j for the regular polygon (right). .
applied to the vertices. Therefore the resulting Hessian matrix H λ is circulant with respect to its 2 × 2 blocks: The spectrum of this block circulant matrix is made of the union of the spectra of the following n matrices of size 2 × 2 where ρ k = exp(ikθ), k = 0, ..., n − 1. Fore more details the reader can refer to [49] and the references therein. One may note that the symmetry of H λ implies that H n−k = H T k . Moreover, the 2 × 2 matrices described in (43) are all Hermitian (and therefore have real eigenvalues).
In the following we assume that the triangulation defining the functions ϕ j in (10) is symmetric and is made of the triangles T j having vertices (0, 0), (cos jθ, sin jθ), (cos(j + 1)θ, sin(j + 1)θ), 0 ≤ j ≤ n − 1. For convenience we may use the notation T + = T 0 , T − = T n−1 (see Figure 3). With these notations it can be seen that for 0 ≤ j ≤ n − 1 we have (44) ∇ϕ j = 1 sin θ Furthermore, in view of the symmetry of the eigenfunction, a simple integration by parts shows that Denote with M 0 , M 1 , ..., M n−1 the blocks of the first line in M λ . Then for ρ k = exp(ikθ) a root of unity of order n we have where R τ = cos τ − sin τ sin τ cos τ denotes the rotation matrix around the origin with the angle τ in the trigonometric sense. By abuse of notation we will use the same notation for the rotation of angle τ around the origin. Recalling the formula (23) we decompose each one of the blocks In the following, we compute separately the matrices B l ρ k = n−1 j=0 ρ j k M l j R jθ , for l = 1, 2, 3. We denote by Id the identity matrix and J = 0 −1 1 0 . The area of P n is |P n | = 0.5n sin θ.
Note that the matrices M 2 j come from the Hessian of the area. Therefore, by straightforward computations we have M 2 Then we have by the symmetry of the eigenfunction that A xx + A yy = λ 1 /n. The fact that the gradients undergo a rotation when transferred from T − to T + implies the matrix equality We find that −A xx sin θ + A yy sin θ + 2A xy cos θ = 0. With the notations above we have It is immediate to see that (∇ϕ 0 · ∇ϕ 1 ) T + = (∇ϕ 0 · ∇ϕ n−1 ) T − = − cos θ|∇ϕ 0 | 2 . Keeping in mind that |P n | = 0.5n sin θ and |∇ϕ 0 | T ± = 1/ sin θ we get Of course, the other blocks on the first line are all equal to zero. Therefore we obtain It can be noted that since A xx + A yy = λ 1 /n and −A xx sin θ + A yy sin θ + 2A xy cos θ = 0 we can deduce that Using these relations and the computations above we find that It remains to compute the contribution of the terms M 1 j . Let us recall that due to the symmetry of the triangulation defining ϕ i we have, denoting U j = (U 1 j , U 2 j ) the solutions of (18) with the normalization (21), that Note that this implies that R T jθ U j = U 0 • R T jθ . For 0 ≤ j ≤ n − 1 we have . Remark 4.8. For 0 ≤ j ≤ n − 1 the sum of the elements which are not on the diagonal of M 1 j R jθ is zero. This is a consequence of the fact that a(U 1 0 , which simply comes from the change of variables y = R T kθ x and the fact that U 2 0 is odd with respect to y and U 1 0 is even with respect to y (see Proposition 4.7).
The next result shows that the eigenvalues of B ρ k , and as a consequence those of M λ , can be expressed in terms of u 1 , U 1 0 , U 2 0 .
Moreover, the eigenvalues of B ρ k are given by As a consequence, the eigenvalues of the Hessian matrix M λ given in (23) are exactly µ j , j = 0, ..., 2n − 1.
Proof. In view of the previous computations we have The formulas follow directly from Proposition 4.7 and Remark 4.8.
Proof: When k = 0 the computations in Proposition 4.10 and the fact that Ω ∇U 1,2 Using the relations computed above we find that α k − γ k = 0. Using the second formula for γ k in Theorem 4.9 and the fact that t y = (0, 1, ..., 0, 1) is an eigenvector of M λ from Proposition 4.6 we find that β k = γ k . The case k = n − 1 follows from B ρ n−1 = B ρ 1 .
Corollary 4.13. We have B ρ k = B ρ n−k (with indices modulo n). Therefore: 1. B ρ k and B ρ n−k have the same eigenvalues.
2. If n is odd then the spectrum of M λ consists of 4 zero eigenvalues and n − 2 double eigenvalues.
3. If n is even then B ρ n/2 is diagonal and the spectrum of M λ consists of 4 zero eigenvalues, n − 4 double eigenvalues and another two eigenvalues that can be found on the diagonal of B ρ n/2 .
For the sake of completeness, in the following we give a short proof that the regular polygon is a critical point for P → |P |λ 1 (P ). This result is known and can be recovered, for instance, using ideas from [19] or [6, Chapter 1]. The proof given below relies on the representation formulas for the gradient given in Theorem 2.2.
Theorem 4.14. The regular polygon is a critical point for x → A(x)λ 1 (x).
Proof. Fix the regular polygon P n inscribed in the unit circle with a 0 = (1, 0) and denote by λ 1 its first eigenvalue. Consider the functions ϕ i , i = 0, n − 1 defined in (10) and suppose they are symmetric like in the right picture in Figure 1. For i ∈ {0, ..., n − 1}, the components 2i, 2i + 1 of the gradient of the objective function are given by In view of the symmetry of the polygon and of the first eigenfunction, it is enough to perform the computations for i = 0. We have ϕ 0 = (1, −1/ tan θ)1 T + + (1, 1/ tan θ)1 T − . This already shows that Using the expression of ∇λ 1 (x) and (45) we find that Using (47) we simplify the above expression to Adding (48) and (49) we find that the first two components of the gradient of x → A(x)λ 1 (x) are zero. By symmetry, all the other components are zero and P n is indeed a critical point.

A priori error estimates for the coefficients of the Hessian matrix
The eigenvalues of M λ are described analytically in the previous section, but the formulae do not allow us to prove that these eigenvalues are non-negative. In view of Proposition 4.5 proving that M λ has 2n − 4 eigenvalues that are strictly positive is enough to infer the local minimality of the regular polygon. In this section we describe how we can certify numerically this fact. In order to achieve this we provide a priori error estimates concerning numerical approximations based on finite elements for α k , β k , γ k given in Theorem 4.9.
First we refer to classical certified estimates for the approximation of the first eigenpair and of the second eigenvalue on the regular polygon P n using P 1 finite elements. In a second step, we get certified estimates for the finite element approximation of the function U i . In the last step we get certified approximation results for the coefficients of the Hessian matrix.

5.1.
Step 1. Certified approximation of the first eigenpair and of the second eigenvalue. In the literature one can find certified approximation for the first eigenvalue in regular polygons (see for instance [31]). We shortly recall of the results of [37,Theorem 4.3].
Let us consider a triangulation T h of P n . In each triangle T i ∈ T h , the ratio between the smallest edge and the middle one L i is denoted α i and the angle between these two edges is τ i . Then, we denote Following [37, Section 2], we introduce the constant where the parameter h dictating the size of the mesh is the size of the median edge. Let us denote V h the finite element space associated to T h with P 1 finite elements. Denote by λ k,h , u k,h the k-th eigenvalue of P n and its associated eigenfunction approximated in V h , solving Results of [38] show that As a direct consequence we have Denoting Π 1,h the Lagrange interpolation operator on the vertices of triangles of T h , for functions u ∈ H 2 (P n ) we have For each u ∈ H 1 0 (P n ) let us denote P h (u) the projection of u onto the finite element space V h , namely the solution of In particular, for u = u 1 ∈ H 2 , using D 2 u L 2 = ∆u L 2 ([22, Theorem 4.3.1.4]), we get In order to estimate the error for the eigenfunction, let u 1,h be an L 2 -normalized, finite element approximation of the first eigenfunction given by (50).
Let us denote by p = P h (u 1 ) and decompose p = αu 1,h + p, where Pn pu 1,h dx = 0, α ∈ R. Note that changing the sign of u 1,h still gives an L 2 -normalized solution, therefore we may assume α > 0 in the previous decomposition.
As we know that Using the Poincaré inequality on the orthogonal of u 1,h in V h , we get We obtain the error estimate We compute the following bounds for p L 2 , ∇p L 2 , which are immediate from the definition of p and the projection operator P h : In order to conclude we need bounds for α. We have Pn p 2 = α 2 + Pn p 2 , which shows that This estimate can be written in a quantitative form using (55) and (54). Since α > 0, for h small enough, an explicit lower bound for α can also be found.
In the same way we obtain the L 2 error estimate for the first eigenfunction It can be noted that the optimal rates of convergence are obtained in (56) and (57). Moreover, the term of order O(h) in (56), which dominates the estimates comes from the interpolation error bound for ∇u 1 − ∇p L 2 while the remaining terms are of higher order O(h 2 ).

5.2.
Step 2. Certified approximation of U j . We begin with some generic approximation results for solutions of the Laplace equation with Dirichlet boundary conditions with singular right hand sides.
From the Gagliardo-Nirenberg interpolation inequality (all norms in R 2 are taken to be the Fourier transform ones), we get Finally, we get the conclusion.
Theorem 5.2. Let U ∈ H 1 0 (P n ) be the solution of Proof. We denote U reg , U sing the solutions of We introduce the auxiliary functions V reg , V sing ∈ V h , V reg = P h (U reg ), V sing = P h (U sing ) the finite element solutions of For V sing , the estimate (58) from Lemma 5.1 holds, and gives ∇U sing − ∇V sing L 2 ≤ f sing while for V reg the estimate from (54) gives Let us denote V = V reg + V sing and defineṼ = V − ( Pn V u 1,h dx)u 1,h . Then we have We have thatṼ is the finite element solution of By the Poincaré inequality in the orthogonal of u 1,h we get Finally, Theorem 5.3. With the notations of Theorem 5.2, the following estimate holds Remark 5.4. It can be seen that the estimates from Theorems 5.2, 5.3 become explicit as soon We present below some inequalities that help obtain upper bounds for all other quantities presented here.
Using the fact that U is orthogonal on the first eigenfunction u 1 we find SinceṼ is the projection of V on the orthogonal of u 1,h in V h we immediately have Ṽ L 2 ≤ V L 2 and ∇Ṽ L 2 ≤ ∇V L 2 . We obtain the following estimates for U h : Since U h is orthogonal on u 1,h the first eigenfunction associated to λ 1,h we have ∇U h This implies Remark 5.5. In practice, the singular right hand side that we consider is of the following type.
Then for every γ ∈ (0, 1 2 ) we have Indeed, for every ϕ ∈ H and use the trace theorem for ϕ from H Practical estimate. In order to estimate g L 2 (S) above, we notice that Here, we have used that ∂u 1 ∂x ∈ H 1 0 (S) for symmetry reasons and angular behavior, together with ∂ 2 u 1 ∂y 2 ≤ 0, from symmetry and convexity of the level lines of u 1 . In order to estimate S u 2 1 dx, we can use the following If the approximation by finite elements has the symmetry of the n-gon, then Analysis of the function U 0 = (U 1 0 , U 2 0 ). As we have seen in Section 4, it is enough to concentrate on the function U 0 defined in (18) with the normalization condition (21). Recall the definition of ϕ i is given in in (10) (see also Figure 1).
Denote by f 0 ∈ H −1 (P n , R 2 ) the right hand side from (18) for j = 0. Since the index 0 is fixed, we shall drop it from U 0 , its right hand side f 0 and ϕ 0 . We denote by T + = T 0 and T − = T n−1 the upper and lower triangles of the support of ϕ. Then ∇ϕ = 1 T + 1, − 1 tan θ + 1 T − 1, 1 tan θ and ∀x ∈ P n , ∇ϕ(x) = 1 sin θ The right hand side f ∈ H −1 (P n , R 2 ) is the distribution given by Recall that S λ 1 = (|∇u 1 | 2 − λ 1 u 2 1 ) Id −2∇u 1 ⊗ ∇u 1 . Using integration by parts, we notice that We observe that in the expression of Pn S λ 1 ∇ϕ the first term cancels for symmetry reasons. We may also use the fact that the regular polygon is critical for λ 1 (P n )|P n | so that Working with a symmetric triangulation for ϕ j (Figure 3) and a mesh that is exact on T j and respects the symmetries of the regular polygon ( Figure 4) the uniqueness of the first discrete eigenfunction u 1,h implies that (46) holds also for the discrete quantities. In particular T + (∂ x u 1,h ) 2 + T + (∂ y u 1,h ) 2 = λ 1,h /n and − sin θ T + (∂ x u 1,h ) 2 + sin θ T + (∂ y u 1,h ) 2 + 2 cos θ T + ∂ x u 1,h ∂ y u 1,h = 0. We denote by and using the previous relations we find that s λ 1,h = −2λ 1,h /n. Therefore |s λ Proposition 5.6. The following inequality occurs Proof: The proof is straight forward by direct computation, taking into account the vector norm inequality (a ⊗ b)c ≤ a b c and the Poincaré inequality We also used the symmetry of the mesh, which gives Practical estimate. We estimate below the quantities needed for the estimates in Theorem 5.2. In order to estimate f i H −1 , f i reg L 2 and f i sing H − 1 2 −γ we use the following notations: Recall that Q(v) is the projection of v on the orthogonal of u 1 . Therefore For the H −1 estimate we work with the formula involving Q(v) and we have Which implies that Let us denote S + , S 0 , S − the segments [0, exp(iθ)], [0, 1], [0, exp(−iθ)] in the complex plane, and n = (n x , n y ) the outside normal of a domain. We have We decompose each term in A i = A i reg + A i sing , the regular part given by the first integral over T + ∪ T − and the singular part given by the sum of the last two integrals over the boundaries of ∂T + and ∂T − . Following [22,Lemmas 3.4 . Therefore, for the regular parts we have We explicit ∇ϕ · n on S 0 , S + , S − and we obtain We have ∂ x u 1 L 2 (S ± ) = cos θ ∂ x u 1 L 2 (S 0 ) and ∂ y u 1 L 2 (S ± ) = sin θ ∂ x u 1 L 2 (S 0 ) since the normal component of the gradient of u 1 is zero on these segments. Therefore, using Remark 5.5 we obtain ∂y + ∂T + v∇ϕ · ∇u 1 n y + ∂T − v∇ϕ · ∇u 1 n y , and we observe that the boundary integrals vanish. For the regular parts we have the same estimates as before It is straightforward to see that C 1 , C 2 are L 2 distributions, C 2 = 0 and C 1 L 2 = |s λ 1 | = 2λ 1 /n. Finally, we get Using the fact that ∂ x u 1 2 T + + ∂ y u 1 2 T + = λ 1 /n we may also use the slightly weaker, but simpler bounds below:

5.3.
Step 3. Estimates for the eigenvalues of M λ . As shown in Theorem 4.9 and Proposition 4.10 the eigenvalues of M λ can be expressed in terms of u 1 and (U 1 0 , U 2 0 ). As we saw in the previous sections, the terms containing derivatives of u 1 can be well approximated using P 1 finite elements using an estimate of order O(h) with explicit constants.
Results of the previous section show that the estimate of the computation error for U behaves like h 1 2 −γ . Trying to bound directly the error for the eigenvalues of M λ will give estimates of the same order, which in practice are not fine enough to provide bounds that allow to certify that the non-zero eigenvalues of M λ are positive.
However, it turns out that the estimate of the coefficients of the shape Hessian matrix of the eigenvalue is better, namely in h 1−2γ , as a consequence of the particular structure of the coefficients. As shown in [20,Section 5] defining and solving an auxiliary problem using the same bilinear form can double the speed of the convergence.
We use the notations of Theorem 5.2 for two generic problems with solutions U a , U b corresponding to the right hand sides f a , f b (not necessarily those explicited in the previous section). As well, we use the associated notations We denote the bilinear forms a : H 1 0 (P n ) × H 1 0 (P n ) → R, a(u, v) = ∇u · ∇v − λ 1 uv, Our objective is to estimate error terms of the type |, in order to get an estimate of order h 1−2γ for α k , β k , γ k in Theorem 4.9. We have We estimate each term of the right hand side, separately, the most delicate being the first one. First term.
As a consequence of Lemma 5.1 applied for the L 2 -norms of both the functions and their gradients we get a control in h 1−2γ .
which, in view of inequality (63), leads to an approximation of order h. Third term.
The last inequality is a consequence of the fact that a h (·, ·) is a scalar product on {u 1,h } ⊥ in V h and of the Cauchy-Schwarz inequality together with the observation that a h (v, v) ≤ |∇v| 2 . Using inequality (64) we get an approximation of order h.
Remark 5.7. The problematic term in the previous estimates can be simplified when the two distributions and associated solutions have opposite parity properties. Indeed, suppose that leading to an estimate of order h 3/2−γ , for γ ∈ (0, 0.5).
Below we show how to choose the functions in the above estimates in order to obtain the desired bounds for the quantities described in Theorem 4.9. Since in the case k = 0 we have α 0 = β 0 = γ 0 = 0 we focus only on the cases 1 ≤ k ≤ n − 1.
Remark 5.8. It can be noted that the error estimates above can already be applied for terms of the type a(U 1,2 j , U 1,2 l ) that appear in the expressions of M λ and α k , β k , γ k . However, if multiple such terms are present in some expression, a direct error estimate will accumulate the errors and the final results will be unusable for reasonably large h. It is best to choose properly the functions U a , U b beforehand and apply the error estimate only once.
Expressing the derivatives of u 1 in the (v j , v j ) basis, by direct computation one gets Consider now the discrete version of f α k , replacing u 1 and λ 1 by their discrete approximations Working under the hypothesis that the mesh T h has the symmetries of the regular polygon and that the triangles T j are meshed exactly, we have (f α k h , u 1,h ) H −1 ,H 1 0 = 0. By direct computation we obtain (cos(j + 1)kθ + cos jkθ) The term β k . Denote by so that a(U 2 0 , W β k ) allows us to express β k (see Proposition 4.10). The orthogonality of U 2 0 on u 1 implies that cos(j + 1)kθ − cos jkθ sin θ T j − cos(2j + 1)θ − sin(2j + 1)θ − sin(2j + 1)θ cos(2j + 1)θ ∇u 1 · ∇v.
Proof: Recall that the estimates given in Section 5.1 allow us to obtain explicit bounds for T 0 (∂ x u 1 ) 2 and T 0 (∂ y u 1 ) 2 of order O(h). Denoting q k = 2n(1 − cos(kθ))/ sin θ we have the following. • We note that (f β k sing , v) H −1 ,H 1 0 = 0 for every function v that is odd with respect to y. Since U 1 0 and its numerical approximation verify this hypothesis as soon as T h is symmetric with respect to the x axis we may apply Remark 5.7 and obtain a better error estimate.
• For γ k = −2|P n |a(U 1 0 , W γ k ,1 ) we apply (66) with U a = U 1 0 , U b = W γ k ,1 . We note that (f γ k ,1 sing , v) H −1 ,H 1 0 = 0 for every function v that is even with respect to y. Since U 1 0 and its numerical approximation verify this hypothesis as soon as T h is symmetric with respect to the x axis we may apply Remark 5.7 and obtain a better error estimate.
• For γ k = 2|P n |a(U 2 0 , W γ k ,2 ) we apply (66) with U a = U 2 0 , U b = W γ k ,2 . In conclusion, the terms α k , β k , γ k admit quantified approximations of order O(h 1−2γ ) for every γ ∈ (0, 1/2). 6. Numerical simulations 6.1. Local minimality. Given the regular polygon P n with n sides inscribed in the unit circle with a vertex at (1, 0), we divide it into n equal slices used in the definition of ϕ i , like in Figure  3. Then we give an integer m ≥ 1 and for each one of the triangles T j , j = 0, ..., n − 1 we construct a mesh T h consisting of congruent triangles similar to 1 m T j . In this way we obtain a mesh with median length h = 1/m. Examples are given in Figure 4. With this definition of mesh T h all triangles in the mesh are similar and the constant C 1 defined in the beginning of Section 5.1 can be explicitly identified in terms of n. Given the mesh T h we compute using P 1 finite elements: • the first two eigenvalues λ 1,h , λ 2,h and the first eigenfunction u 1,h of the discrete Dirichlet-Laplace eigenproblem (50).
for the discrete distributions f j h , j = 0, ..., n − 1 given by using the normalization . Therefore we obtain the approximations of α k , β k , γ k from Theorem 4.9 that are of order O(h 1−2γ ) for γ ∈ (0, 1), with explicit error bounds given in the previous section. The procedure described above provides for each k = 1, ..., n − 1 intervals I α k , I β k , I γ k for which we have the guarantee that α k ∈ I α k , β k ∈ I β k , γ k ∈ I γ k . Using the interval arithmetic toolbox Intlab [44] we find intervals I j (h, γ) containing the eigenvalues µ j , 0 ≤ j ≤ 2n − 1 of M λ described in Theorem 4.9. Given a value of h and the associated numerical approximations we obtain a whole range of intervals I j (h, γ) for γ ∈ (0, 0.5). Note that changing γ at fixed h is not a difficulty since this parameter appears only in the choice of constants and exponents. When γ is close to zero we obtain a weak estimate in Theorem 5.2 while for γ close to 0.5 the constants in the estimates from Remark 5.5 become very large. An appropriate choice for γ is made using a simple grid search. If among the intervals I j (h, γ) we obtain only two that contain zero then we conclude, based on Proposition 4.5, that the regular polygon P n is a local minimum for P → |P |λ 1 (P ). If this is not the case we decrease h and we repeat the procedure. Remark 6.1. Numerical algorithms employed in scientific computing use floating point arithmetic. As a consequence there is a difference between the exact discrete solution of the finite element problem and the one given by the numerical algorithm. The sources of error are as follows: • the numerical mesh is a slight perturbation of the exact mesh, leading to perturbations in the mass and rigidity matrices. • the linear systems are solved using iterative methods with a stopping criterion related to the residual vector. In general, it is admitted that errors coming from the above considerations are smaller than the theoretical error estimates shown in Theorem 5.10. The condition number of the linear systems involved is of order O(h −2 ), therefore, we expect that for h ≥ 10 −4 the machine errors do not dominate in the estimation of λ 1 . Moreover, for the gradient terms and for U 0 , which have an  Table 1. Size of the computational problems for the finite element computations.
even weaker convergence rate, the discretization error is expected to dominate machine errors. We also make this assumption in the following.
Formulas given in Theorem 4.9 and Proposition 4.10 allow us to compute the eigenvalues of the Hessian matrix in knowing the first eigenfunction u 1 on P n and the pair (U 1 0 , U 2 0 ) solution of (18). Using P 1 finite elements it is straightforward to approximate the first eigenpair. Given a mesh T h with N v vertices, and denoting by (φ i ) Nv i=1 the P 1 basis functions, the rigidity and mass matrices are defined by The first eigenpair and the second eigenvalue are approximated by solving the generalized eigenvalue problem Ax = λBx. Denote by x 1 the eigenvector associated to the first eigenvalue. Then (61) is solved by considering embedding the orthogonality on u 1,h in the linear system: The constraint vector c is given by c = x T 1 M and the right hand side f is computed by evaluating (f 1,2 0 , φ i ) H −1 ,H 1 for every φ i in the finite element basis. In order to have an error estimate small enough such that the interval around the eigenvalue does not contain zero rather small values of h need to be considered, leading to large computational problems. The value of h and the number of degrees of freedom (d.o.f) for the computational problems are listed in Table 1. Therefore, in order to be able to solve these problems the software FreeFEM [23] is used in its parallel version together with the libraries PETSc [4], SLEPc [28], Hypre [17]. The computations use 200 processors and are run on the cluster Cholesky from the IDCS Mesocenter at Ecole Polytechnique. The error estimates allow us to obtain sufficiently small intervals for h = 10 −4 for n ∈ {5, 6, 7, 8}. The resulting eigenvalues and quantities needed are given to the interval arithmetic library Intlab [44]. The library is then used to compute the interval enclosures for the eigenvalues. The non-zero eigenvalues and the corresponding enclosures are given in Table 2. The results shown in Table 2 indicate that the regular polygon is a local minimizer for problem (2) for n ∈ {5, 6, 7, 8}. In Table 3 we estimate the largest mesh size h for which the certified numerical computations validate the local minimality of the corresponding regular polygon. Exploiting the symmetry of the eigenfunction and of the functions U 1 0 , U 2 0 the size of the problems can be further reduced in half. Remark 6.2. The results shown in this section prove the local minimality of the regular polygon when neglecting errors coming from floating point computations. Most algorithms are designed such that these errors are minimized and therefore it is generally agreed that these errors are smaller than the errors between the continuous solution and the exact discrete one. However, guaranteeing that the floating point errors are small enough it is a non-trivial matter that needs to be addressed in future works. Ideally, the whole computation of the finite element problems should be handled using an interval arithmetic library like Intlab [44], which is a non-trivial task in view of the minimal size of the problems listed in Table 3.
It is possible to compute the eigenvalues of the Hessian matrix for higher n, without guarantee that the numerical eigenvalues are precise enough. Nevertheless, it is well established that a priori estimates are rather pessimistic and the following results might precise enough. In Table  4 we present the non-zero eigenvalues of the Hessian matrix for h = 10 −3 for 9 ≤ n ≤ 15. These  Table 3. Approximately optimal mesh sizes and number of degrees of freedom for which currently known a priori estimates allow to certify the local minimality. eigenvalues are positive, suggesting that the regular polygon is still a local minimzier in these cases.
6.2. General gradient descent simulations. The gradient of the first eigenvalue with respect to the coordinates of the vertices is given in Theorem 2.2. Using these formulas is straightforward to implement a gradient descent algorithm starting from random initial polygons.
Simulations were preformed for the minimization of the first eigenvalue for n ∈ [5,15] and in every case the result of the optimization was a polygon very close to being regular. In order to see how close to being regular is the polygon ω n given by the simulation the following information is given in Table 5: the optimal numerical first eigenvalue, the difference between the maximal and minimal edge lengths, the difference between the maximal and minimal angles (in radians), the difference between the optimal numerical eigenvalue and the precise first eigenvalue of the regular polygons ω * n given on the following web page: http://hbelabs.com/regularpolygon/ index.html (based on the article [32]). Repeating the simulation starting from random initial polygon always gives similar results. The results shown in Table 5 indicate that the optimal numerical polygons ω n found by the numerical algorithm are close to being regular. Furthermore, the value of the objective function is as close to the precise value given for the actual regular polygon ω * n , as the precision of the numerical computations allows. These computations further suggest that the regular polygon is indeed the global minimizer for (2). 3.1e-3 1.7e-6 Table 5. Results of the gradient descent optimization algorithm with random initial polygons.

Reduction of the proof of the conjecture to a finite number of numerical computations
In this section we provide a strategy for proving the conjecture using a finite number of computations for a given number of sides. This strategy works under the implicit assumption that the conjecture is true! In order to justify that for every n the conjecture can be reduced to a finite number of numerical computations, we begin with some theoretical analysis. Assuming the area of a polygon with n sides is fixed (say π), we shall find a value D max such that if the diameter of the polygon exceeds D max then the polygon cannot be optimal for (2). As well, we shall find a minimal value for the length of the edges e min and for the inradius r min of an optimal polygon. All these results (which depend on n), produce a compact set of polygons (seen as subset of R 2n−4 ) outside which any polygon cannot be optimal for |P |λ 1 (P ).
We denote by P n the closure of the class of simple polygons with at most n edges for the Hausdorff distance of the complements. A polygon belonging to this class may be degenerate in the sense that one vertex can belong to a different edge. Depending on how this occurs, this may lead to a disconnection, i.e. a union of two polygons. However, as soon as a polygon is optimal, disconnection can not occur.
Let us denote for every n ≥ 3 the minimal value for the scale invariant formulation by l * n = min{|P |λ 1 (P ) : P ∈ P n }. It is known that l * n < l * n−1 (see [24,Section 3.3]).
Theorem 7.1. Let n ≥ 3. There exists a value D max > 0 such that if P ∈ P n , |P | = π and diam(P ) > D max then πλ 1 (P ) > l * n . In other words, when searching the minimizer in the class of n-gons of area π, it is enough to restrict to polygons with diameter less than or equal to D max . This information is crucial in order to limit the number of numerical computation and leads to a formal, inductive, proof of the conjecture. The value of D max can be computed and depends on l * n−1 and λ 1 (P n ). Proof: The proof is inspired by the surgery argument of [11], where the authors propose a precise way to estimate the diameter of an optimal set in relationship with the first eigenvalue.
The key idea is that if the diameter of an optimal set is too large, one can cut the set with a strip of positive width in order to produce a better one. The main difficulty in our case is that cutting a polygon having n edges with a strip may produce a union of polygons, some of which may potentially have more than n edges, making them non-admissible. In order to handle this situation, further analysis is necessary.
Then, the solution of this problem is the same as the solution of the constrained problem with area π set in (2). Let us denote by Q n an optimal polygon, having area π. Let K ≥ l * n /π be fixed. For instance, K may be obtained using a numerical approximation from above of λ 1 (P n ).
Surgery. In order to get the bound on the diameter, we shall use the surgery results of [11]. Let us set the following constant in the class P n . We recall that the torsion energy of P is defined by E(P ) = min Indeed, if for some P ∈ P n , P ⊆ Q n we have E(P ) + c|P | < E(Q n ) + c|Q n |, then from [11, Lemma 3.1]) we would get |P |λ 1 (P ) < |Q n |λ 1 (Q n ), in contradiction to the optimality of Q n .
Step 2. (Use of [11,Corollary 4.3]) Let w be the torsion function of Q n . Assume a ∈ R and denote by S r (a) = {(x, y) : r − a < x < r + a} an open strip in R 2 . Assume that the interior of the strip intersects Q n and does not contain any vertex. In this case, the intersection of the strip with Q n is a union of trapezes {T j } j∈J . When removing any of these trapezes, one splits the polygon Q n in two (or more, if a vertex is on the boundary of the strip) polygons.
In fact, taking a closer look to the argument of [11,Corollary 4.3], leads as well to As a consequence of [11, Lemma 3.1] this implies This last inequality leads to a contradiction of the optimality of Q n only if the open set Q n \ T j consists in a union of polygons, each one with at most n edges. In this case, it is enough to pick the one with minimal first eigenvalue and contradict the optimality of Q n . Of course, it may happen that one of the connected components of Q n \ T j is a polygon with more than n edges, as new edges could be produced by the surgery procedure. We shall prove that if the diameter is larger than some computable constant, then there exists some suitable strip S r (a) and a suitable trapeze T j such that each connected component of Q n \ T j is a polygon with at most n edges. This contradicts the optimality of Q n .
Step 3. (Preparatory facts) We know from the Saint-Venant inequality that The following results is, for instance, contained in [11, Lemma 2.2]: if w(x 0 ) ≥ η > 0, then Consequently, if we consider a strip S 2r 0 (a) such that max S 2r 0 (a) w > C 2 0 , then, recalling that C 0 = r 0 , where x 0 is a maximum point of w in S 2r 0 (a). In particular Let us introduce the natural number ( · denotes the integer part) Clearly, if the diameter of Q n is larger than 8C 0 k, then taking the x-axis along the diameter, there will be at least one strip of width 4C 0 where the mass of w is less than C 2 0 .
We recall now the following inequality, for which we refer to [50]. Let Ω be a bounded, open simply connected set in R 2 . Let w Ω be the torsion function in Ω. We have w Ω (x) = Ω G Ω (x, y)dy, where G Ω (x, y) is the Green function for the Dirichlet-Laplace operator on Ω. From the Cauchy-Schwarz inequality we have In [50, Proof of Theorem 1.5, inequality (5. 16)] it is shown that if πR 2 0 = |Ω| then where d(x) is the distance from x to ∂Ω. This leads to the estimate We use this inequality for Ω = Q n , so that |Q n | = π, getting the bound w Qn (x) ≤ 2 √ 2d(x) 1/2 . We introduce now e * , d * such that and d * = π e * . Lemma 7.2. The diameter of Q n can not be larger than 2d * + (k + n − 2)8C 0 .
Proof. Assume for contradiction that there are two vertices a 0 , a m such that the diameter of Q n is the segment [a 0 , a m ] and that its length is larger than 2d * + (k + n − 2)8C 0 . Around the midpoint of [a 0 a m ] we build k + n − 2 adjacent strips of width 8C 0 . Outside the strips there are two sub-segments of [a 0 a m ], each having length at least d * (see Figure 5). We remove at most n − 2 strips having a vertex in their interior and among the remaining k strips there is one, say S 4C 0 (a) such that max From the choice of the strip, the set S C 0 (a) does not contain any vertex of the polygon Q n , so that an edge either crosses the strip from one side to the other, or it stays on the same side. In particular, this implies that Q n ∩ S C 0 (a) is a union of open trapezes {T j } j . Moreover, we get for each such trapeze |Q n \ T j |λ 1 (Q n \ T j ) < |Q n |λ 1 (Q n ).
Assume we remove one trapeze, say T j , from Q n ∩ S C 0 (a) and get two polygons P l T j and P r T j , which together have n + 4 edges. There are two possibilities.
(1) Both polygons P l T j and P r T j have no more than n edges. This situation contradicts the optimality of Q n .
(2) One of P l T j and P r T j has n + 1 edges and the other one has 3 edges. In the following, we suppose that the second situation above occurs for each trapeze T j , otherwise we contradict optimality. We claim that on one side of the strip there are only triangles.
If there is only one trapeze, there is nothing to prove. Assume for contradiction that there are two trapezes, which when removed generate triangles on both sides of the strip. From simple connectedness, there is a continuous curve contained inside the polygon, joining the interiors of the two triangles. See Figure 5 (a). This curve crosses the strip at least one more time, implying the presence of at least another trapeze, which cannot leave a triangle on either side when removed without disconnecting the polygon. Therefore, removing this trapeze, we split the polygon in two polygons with less than n edges contradicting optimality.
In conclusion, removing any one of the trapezes T j generates triangles, all situated on one side of the strip. Assume this occurs on the left. Now, we choose the triangle containing the vertex a 0 on the left, which is at distance at least d * from the strip and the trapeze which isolates it in a triangle. We continuously move the strip S C 0 (a) to the right (and the trapeze with it) up to the moment when the strip touches a first vertex. This vertex can be a neighbor of a 0 ( Figure  S  5 (a)) or a different vertex a k (Figure 5 (c)). In any case, the trapeze will split the polygon in either two or three polygons and the number of edges for each polygon is at most n. Moreover, one polygon is the triangle with a vertex in a 0 . The area of this triangle together with the trapeze is at least d * e/2, where e is the length of the longest vertical edge of the trapeze, on the right side of the strip. This set is fully contained in the polygon, so has area at most π, meaning that e 2 ≤ e * . Using inequalities (71)-(72) we get that the maximum of w Qn on the trapeze is below C 2 0 . This contradicts the optimality of the polygon. Theorem 7.3. Assume P = [a 0 ...a n−1 ] ∈ P n is such that |P | = π, diam(P ) ≤ D max . There In other words, an optimal polygon of area π in P n can not have an edge smaller than a certain threshold. To observe this fact, it is enough to choose δ 0 such that (74) l * n−1 − Cδ 1 2 0 > l * n . The following type of result has been proved by Davies in [16] and refined by Pang in [41]. We give a short proof below, based on the comparison with the torsion function.
Lemma 7.4. Assume P = [a 0 ...a n−1 ] ∈ P n is such that |P | = π, diam(P ) ≤ D max . Let Q ∈ P n , Q = [b 0 ...b n−1 ] such that for every i = 0, . . . , n − 1, |a i b i | ≤ δ. Then Proof. Assume in a first step that Q ⊆ P and denote w Q , w P the associated torsion functions. The inequality is a consequence of [11, Inequality (2.6)] which gives and of the estimate which is a consequence of (71) applied to w P and of the harmonicity of w P − w Q on Q.
In general, if Q ⊆ P , we use the previous argument and compare both λ 1 (P ), λ 1 (Q) with λ 1 (P ∩ Q).
Assume at least one of the angles a 0 , a 1 is convex, for example a 0 . Then we move the point a 0 towards a 1 continuously, denoting it a t 0 = (1 − t)a 0 + ta 1 . If the segment [a n−1 a t 0 ] does not meet any other vertex of the polygon for any t ∈ (0, 1), then we denote Q the new polygon obtained for t = 1. Clearly, Q ∈ P n−1 and Lemma 7.4 can be applied to get From Makai's inequality [39] we know that λ 1 (P ) ≥ 1 4ρ 2 P , where ρ P is the inradius. Since On the other hand, ρ Q ≥ ρ P − δ, hence Finally we get and we conclude this case. If for some t ∈ (0, 1) the edge [a n−1 a t 0 ] meets a vertex. Then the inequality above is still true, and the polygon P t is split in two polygons, each one with at most n − 1 edges. See Figure 6 (left). We choose the one which has the lowest eigenvalue and repeat the previous argument.
If both angles a 0 , a 1 are concave, we consider the same type of movement a t 0 = (1 − t)a 0 + ta 1 . If the segment [a n−1 a t 0 ] does not meet any other vertex of the polygon for any t ∈ (0, 1), then we denote Q the new polygon obtained for t = 1 and follow the previous argument. The difference occurs if [a n−1 a t 0 ] meets a vertex, say a k . In this case the angle a k is convex. However, there is no splitting in this case. We continue the movement, moving at the same time a k parallel to [a n−1 a 0 ], denoted a s k See Figure 6 (right). If the movement finishes at t = 1, we apply previous argument. If the movement blocks because the segments [a k−1 a s k ], [a s k , a k+1 ] touch another vertex then the polygon splits, since we move a convex angle towards the interior of the polygon, and we stop following the same argument as in the previous case. If it blocks because [a n−1 a t 0 ] meets another vertex, we treat it the same way as a k and continue the movement. The blocking can occur at most n − 3 times.
Below we show that a strategy to prove the conjecture by a finite number of numerical computations can be followed, provided of course that the conjecture is true. However, from a practical point of view, this strategy is far from being optimal.
Theorem 7.5. Provided the conjecture is true, for every n ≥ 5 its proof can be reduced to a finite number of numerical computations.
Proof. We know the inequality is true for n = 4. Assume now that the inequality is true for polygons with up to n − 1 edges. Then recalling the notation P n = [a * 0 ...a * n−1 ] for the regular polygon inscribed in the unit circle with one vertex at a * 0 = (1, 0) we have l * n−1 = λ 1 (P n−1 )|P n−1 |. We have a certified estimate from above and from below for this value. In order to prove the conjecture for n edges, we shall analyze problem (68) in the following steps.
Step 1. Compute a certified approximation of the first eigenpair (λ 1 , u 1 ) on P n . The certified approximation of the eigenfunction u 1 holds in H 1 0 (P n ).
Step 2. For the regular polygon P n inscribed in the unit circle, having the vertex a 0 = (1, 0) we compute the spectrum of the shape Hessian of λ 1 (P n )|P n |, and we certify the positivity for 2n − 4 of its eigenvalues using results from Sections 4, 5.3. This concludes that the regular polygon is a local minimum.
From now on, we identify polygons P = [a 0 a 1 . . . a n−1 ], with a i = (x i , y i ), by a point in R 2n−4 having coordinates (x 2 , y 2 , . . . , x n−1 , y n−1 ). We consider the first two points fixed: a 0 = a * 0 , a 1 = a * 1 . Without restricting generality, we can assume that the edge [a 0 a 1 ] is the longest edge in the polygon P . Let us denote P the family of such polygons of P n , identified as a compact subset in R 2n−4 .
Precisely, for a value ε 0 > 0 we have |P |λ 1 (P ) ≥ |P n |λ 1 (P n ) for every P = [a 0 . . . a n−1 ], with a 0 = a * 0 , a 1 = a * 1 , such that ∀i = 2, n − 1, |a i a * i | ≤ ε 0 . Of course, in order to obtain ε 0 , the availability of the constants C and ϑ in Theorem 3.18 is assumed. Let us denote L n this neighbourhood, which is a closed set.
Step 4. Using Theorem 7.1 find an estimate for the minimal measure of an optimal polygon in the class P. Here we use the fact that the maximal length of an edge is precisely [a 0 a 1 ]. Then, we get |P n | |P | |a * 0 a * 1 | ≤ D max .
Using Makai's inequality we get a lower bound for the inradius of an optimal n-gon (called ρ min in the sequel), since 1 ρ 2 P |P | ≤ |P |λ 1 (P ).
In particular, if ρ 2 P ≥ |P |/l * n then P is cannot be optimal. Using Theorem 7.3, we obtain a lower bound on the shortest edge, e min .
All these three geometric constraints: measure, inradius and shortest edge generate a smaller compact set P , defined by purely geometric constraints, such that P ⊆ P in R 2n−4 . In particular, the lower bound on the inradii of such polygons, makes that the inequality in Lemma 7.4 becomes uniform. Below we work with δ ≤ ρ min 4 < 1. Indeed, for every P, Q in the class P the value λ 1 (B ρ min −δ ) is an upper bound for λ 1 (P ), λ 1 (Q), λ 1 (P ∩ Q): for P, Q ∈ P , there exists a universal constant K (with explicit value, issued from Lemma 7.4) such that if the distance between the respective vertices is at most δ then |λ 1 (P ) − λ 1 (Q)| ≤ Kδ The variation of the area ||P | − |Q|| is also controlled by a term of the form K δ, with K = nD max + nπ. There exist universal upper bounds for the first eigenvalue and for the area in P . Therefore, there exists an explicit constant K such that if the distance between the respective vertices of P, Q ∈ P is at most δ then (75) |λ 1 (P )|P | − λ 1 (Q)|Q|| ≤ |P |(λ 1 (P ) − λ 1 (Q)) + λ 1 (Q)||P | − |Q|| ≤ K δ 1 2 .
Step 5. Suppose the conjecture is true and P n is the only minimizer for P → |P |λ 1 (P ) in P . Then, in view of Step 2. above, where a local minimality neighborhood was identified around P n , there exists ε 1 > 0 such that outside the neigbourhood L n of P n in R 2n−4 we have |P |λ 1 (P ) > l * n + ε 1 . If such an ε 1 cannot be found, then a minimizing sequence which does not converge to P n could be constructed, contradicting the hypothesis that the conjecture is valid.
Consider δ > 0 such that K (2δ) 1 2 < ε 1 /4, with K from (75). Moreover, suppose that 2δ ≤ δ 0 with δ 0 from (74). We cover the compact set P \ L n with at most c 2n−4 Choose one of the balls B j enumerated above. Take an admissible polygon P ∈ P \ L n having coordinates (x 2 , y 2 , ..., x n−1 , y n−1 ) in the ball B j . If such a polygon does not exist, there is nothing to be done and we move to the next ball. We evaluate |P |λ 1 (P ) numerically, obtaining a certified estimate interval of length at most ε 1 /4. If this certified computation gives (76) |P |λ 1 (P ) ≥ l * n + ε 1 2 , then P is not optimal and, in view of the choice of the constant δ, no other optimal polygon exists having coordinates in the same ball. If (76) holds for every ball B j containing an admissible polygon then the conjecture is solved. However, the value of ε 1 is not known. For this reason, we start with a value ε 1 = 1 and perform the computations enumerated above. If inequality (76) holds every time there is an admissible polygon in one of the balls then the proof succeeded and we stop. If for some polygon the inequality fails, we divide ε 1 by 2 and restart the computation, etc. This procedure stops in a finite number of steps. Note that from practical point of view, this procedure is completely inefficient, but formally leads to the conclusion. The Saint-Venant inequality states that the maximum of the torsional rigidity among all sets of area π is achieved on the disc. Pólya and Szegö have also conjectured in 1951 (see [42, page 158]) the following.
Conjecture. The unique solution to problem (78) is the regular polygon with n sides and area π.
All the results we have obtained for the eigenvalue transfer similarly to the conjecture above. However, this conjecture is computationally less challenging than the eigenvalue. In particular, there is no additional normalization and orthogonality constraints for w and for the associated material derivatives. The proof of the local maximality goes through the computation of the Hessian matrix of (77) on the regular polygon. The expression of its coefficients was obtained by Laurain in [36]. Recalling that the functions ϕ i are constructed in (10), one introduces the functions U i ∈ H 1 0 (P, R 2 ), i = 0, ..., n − 1 (79) P DU i ∇v = P −(∇ϕ i ⊗ ∇w)∇v + 2(∇w ∇v)∇ϕ i + P v∇ϕ i , for every v ∈ H 1 0 (P ).